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C/3 . ABSTRACT 

An increase in the central density of a neutron star may trigger a phase transition 
from hadronic matter to deconfmcd quark matter in the core, causing it to collapse 
. to a more compact hybrid-star config uration. We present a study of this, building 

■ on previous work by (Lin et all (|2006h . We follow them in considering a supersonic 

| phase transition and using a simplified equation of state, but our calculations are 

. general relativistic (using 2D simulations in the conformally flat approximation) as 

C ' compared with their 3D Newtonian treatment. We also improved the treatment of 

the initial phase transformation, avoiding the introduction of artificial convection. As 
\Q ■ before, we find that the emitted gravitational-wave spectrum is dominated by the 

' fundamental quasi-radial and quadrupolar pulsation modes but the strain amplitudes 

OO , are much smaller than suggested previously, which is disappointing for the detection 

prospects. However, we see significantly smaller damping and observe a nonlinear 
| mode resonance which substantially enhances the emission in some cases. We explain 

the damping mechanisms operating, giving a different view from the previous work. 
Finally, we discuss the detectability of the gravitational waves, showing that the signal- 
to-noise ratio for current or second generation interferometers could be high enough 
to detect such events in our Galaxy, although third generation detectors would be 
needed to observe them out to the Virgo cluster, which would be necessary for having 
a reasonable event rate. 

Key words: hydrodynamics - relativity - methods: numerical - stars: neutron - 
stars: pulsations - stars: phase-transition - stars: rotation 
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1 INTRODUCTION 

The existence of compact stellar objects partially or totally 
consisting of matter in a deconfmed quark phase was al- 
ready predicted long ago (|Bodmerl l 19711 : lltohfll97d : Iwittenl 
1 19841 ) . Such stars are thought likely to originate as a result 
of the conversion of purely hadronic matter in the interior 
of a neutron star (NS) into a deconfined quark matter phase 
when th e density exceeds a certain thre shold (for a review 
see, e.g. IWeberlll999l : lGlendenningll2002l ). We focus here on 
cases where the conversion occurs only in the core of the 



NS, while the outer parts remain unchanged. The theory 
of dense nuclear matter predicts that such compact stars 
[hereafter we will use the term "hybrid quark star" (HQS) 
when referring to these objects] would generally be more 
compact than the progenitor standard NS, and their equi- 
librium radii would be smaller by up to 20%. The potential 
energy W released by the phase transition is expected to be 
of order 
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where R and M are the typical radius and mass of the NS, 
respectively, and AR is the decrease in radius. 
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A first-order phase transition is expected to be the 
most interesting case as far as the dynamics and structure 
of the star are concerned. Such a transition would proceed 
by the conversion of initially metastable hadr onic matter in 
the core int o the new deconfined quark phase l|Zdunik et al.l 
120071 1 2008). The metastable phase could be formed as the 
central density of the NS increases due to mass accretion, 
spin-down or cooling. This could happen soon after the 
birth of the NS in a supernova or it could occur for an 
old NS accreting from a binary companion. The first of 
these channels is likely to give the higher event rate but the 
second is also interesting. In the widely accepted scenario 
for the formation of millisecond pulsars, where an old NS 
is spun up by accretion from a binary companion, the 
amount of spin-u p would be directly r elated to the amount 
of mass accreted (|Burderi et al.|[l999l ). meaning that there 
would be a population of rapidly rotating NSs with rather 
high mass (and hence high central densities). Recent 
observational data seems to confirm that a significant 
proportion of milli second pulsars do indeed tend to be 
high-mass objects (jFreire et al.l 120081 ) which would then 
be candidates for undergoing (or having undergone) a 
phase-transition-induced collapse of the type which we are 
discussing here. While some significant proportion of the 
potential energy release given by Eq. ([1} would probably 
go into neutrino emission, a significant proportion might 
also go into pulsations of the newly-formed HQS, if the 
conversion process to the new phase is rapid enough, 
and this could b e an interesting source o f gravitational 
wave s (GWs; se e iMarranghello et alJ 120021 : iMiniutti et al.l 
120031 ; iLin et all l2006f ). If detected, these GW signals, 
and in particular the identification of quasi-normal mode 
frequencies in their spectrum, could help to constrain the 
properties of matter at the high densities encountered 
here. For non-rotating cold NSs with various composi- 
tions, the related theory of astero seismology has already 
been form ulated in recent years jAndersson fc Kokkotasl 



1998; 



Kokkotas. Apostol atos fc Anderssonl l200ll : 



Benhar. Ferrari fc Gualtieril l2004) 



A problem for making any detailed studies of phase- 
transition-induced collapse is that the description of the 
physics of the phase transition remains very uncertai n 
and controv e rsial (for a recent review see lHorvathll2007l) . 
iDrago et all (|2007t) discussed possible modes of burning of 
hadronic matter into quark matter in the framework of rel- 
ativistic hydrodynamics using a microphysical equation of 
state (EoS). They found that the conversion process always 
corresponds to a deflagration and never to a detonation. 
They also argued that hydrodynamical instabilities can de- 
velop at the burning front. They estimated the correspond- 
ing increase in the propagation velocity of the phase tran- 
sition and noticed that, although the increase is significant, 
it is not sufficient to transform the deflagration into a det- 
onati on in essentially all realistic scenarios. On the other 
hand, iBhattacharvva et all (|2006f ) considered the transition 
as being a two-step process, in which the hadronic matter is 
first converted to two-flavour u and d quark matter, which 
then subsequently transforms into strange quark matter (w, 
d and s quarks) in a second step. They used the relativistic 
hydrodynamic equations to calculate the propagation veloc- 
ity of the first front and found that, in this first stage, a 
detonation wave develops in the hadronic matter. After this 



front passes through, leaving behind two-flavour matter, a 
second front is generated which transforms the two-flavour 
matter to three-flavour matter via weak interactions. The 
timescale for the second conversion is ~ 100 s, while that 
for the first step is about 1 ms. 

Against this background, ILin et aL (2006: hereafter re- 
ferred to as LCCS) carried out a first study of GW emission 
by a phase-transition-induced collapse of a rotating NS to a 
HQS, using a very simplified model for the microphysics and 
treating the phase transition as occurring instantaneously. 
On the basis of 3D calculations using Newtonian gravity 
and hydrodynamics, they obtained waveforms of the emit- 
ted GWs for several collapse models, and found that the 
typical predicted dimensionless GW strain h ranged from 3 
tol5xl0" 23 for a source at a distance of 10 Mpc. The corre- 
sponding energy E gv/ carried away by the GWs was found to 
be between 0.3 and 2.8 x 10 51 ergs. They also determined the 
modes of stellar pulsation excited by the collapse and showed 
that the spectrum of the emitted GWs was dominated by the 
fundamental quasi-radial and quadrupolar pulsation modes 
of the final star. They suggested that the damping of the 
stellar pulsations observed in their calculations was due to 
the production of shock waves leading to the development of 
differential rotation, which proceeds on a timescale of about 
5 ms for typical collapse models. 

The study by LCCS treated the conversion process as 
being instantaneous, and this was mimicked by replacing, 
at the initial time, the EoS describing the hadronic nuclear 
matter in the core (with which the initial equilibrium model 
had been built) by one describing a central zone of decon- 
fined quarks surrounded by a region of mixed phase. The 
material outward of this remained in the original hadronic 
phase. In a first-order phase transition that proceeds via a 
detonation, the conversion front propagates supersonically 
with respect to the matter ahead of it. The sound speed in 
the stellar interior typically has a value of the order of 0.3 
to 0.5 c. Assuming that the quark-matter core has a radius 
of about R ~ 5 km, a rough estimate for the timescale of 
a supersonic conversion gives r ~ R/v ~ 0.05 ms. Clearly, 
this value is not much smaller than the dynamical timescale 
of a NS and so one does not know, a priori, the impact that 
the properties of the conversion would have for the subse- 
quent stellar dynamics taking place on similar timescales. 
Therefore, one should first check how the dynamics of the 
star after the phase transition would depend on the finite 
propagation velocity of the front, and we do this here. 

LCCS represented the hadronic matter by means of a 
polytropic EoS, initially having an adiabatic index 7 = 2. 
When the phase transition was triggered (by changing the 
EoS in the central regions, reducing the pressure support), 
they also replaced the original hadronic EoS by a softer 
ideal-gas type of EoS (with 7 then ranging from 1.95 to 
1.75 depending on the model), which artificially lowered the 
pressure outside the deconfined quark matter core as well 
[see their Eq. (44)]. The lowering of the adiabatic index also 
in the outer regions leads to an increased release of gravi- 
tational binding energy and is hard to motivate on physical 
grounds. 

Although the study by LCCS was an important step 
forward, it was clearly using an extremely simplified model 
for the physics of the NS matter and of the phase transi- 
tion (as well as not including emission processes other than 
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GW emission) and one must question how closely it rep- 
resents such processes occurring in real NSs. Since the full 
problem is complex and involves some input physics which 
is very incompletely known at present, the only way to make 
progress at this stage is, indeed, to work in terms of simpli- 
fied models but making improvements where possible. The 
strategy of this paper is to take the work of LCCS as a 
starting point and then to take some steps forward in the 
direction of including further aspects. In particular, the ef- 
fects of general relativity (GR) certainly play an important 
role in such studies: for instance, the total rest mass of a 
typical Newtonian NS may be almost twice as large as that 
of its counterpart in GR with the same central density and 
the same EoS. The impact of GR is even more pronounced if 
the star is near to the maximum mass limit. Therefore, one 
expects the Newtonian and GR descriptions of NS collapse 
to differ significantly, thus necessitating a proper treatment 
of the GR effects. 

In this study we extend the previous work of LCCS in 
a number of ways: (i) we take into account the effects of 
GR; (ii) we modify the EoS of the stellar matter only in 
regions which undergo the phase transition and take care 
to avoid introducing spurious convection; (iii) we consider a 
larger set of physical models for different values of the stel- 
lar parameters and different properties of the stellar mat- 
ter; (iv) we study the effect of introducing a finite timescale 
for the initial phase transformation which destabilises the 
core rather than treating it as occurring instantaneously. 
Our model for the HQS is based on that of LCCS. The nu- 
merical hydrodynamics simulations are here performed in 
axisymmetry using the conformal flatness approximation to 
GR. While our model remains extremely simplified, we be- 
lieve that these modifications do represent valid and signifi- 
cant steps forward. We feel, however, that we should caution 
against making further elaborate extensions of the hydro- 
dynamics (e.g. to 3D GR hydrodynamics or GR magneto- 
hydrodynamics) until such time as one has the possibility of 
including a greatly improved treatment of the microphysics 
and emission processes so that the models come closer to 
those of real NSs. 

This article is organized as follows: in Section[2]we sum- 
marize the numerical methods used, while in Section [3] we 
introduce the models investigated. In Section U we validate 
our code by performing a comparison with the Newtonian 
models of LCCS. In Section we discuss the results of our 
simulations of GR models and their dependence on the var- 
ious parameters, and in Section [6] we conclude with a sum- 
mary. In Appendix [X] we examine the impact on the col- 
lapse dynamics of considering a finite timescale for the ini- 
tial phase transformation, and in Appendix|B]we present our 
method for determining the damping times for the emitted 
gravitational radiation waveforms. 

Unless otherwise noted, we choose geometrical units for 
all physical quantities by setting the speed of light and the 
gravitational constant to one, c = G = 1. Latin indices run 
from 1 to 3, Greek indices from to 3. 



bv lKomatsu. Eriguchi fc Hachisul dl989aT) (KEH hereafter) 
as im plemented in the code RNS ( Stergioulas fc Friedman! 
Il995l ). This code solves the GR hydro-stationary equations 
for rotating matter distributions whose pressure obeys an 
EoS given by a polytropic relation [see Eq. (fTTj) below] . The 
resulting equilibrium models, which we choose to be rotating 
uniformly, are taken as initial data for the evolution code. 

The time-dependent numerical simulations were 
performed with the CoC o NuT c ode developed by 
iDimmelmeier. Font fc Mulieil (|2002al lbh with a met- 
ric solver based on spe ctral methods as described in 
IDimmelmeier et al.l (|2005l ). The code solves the GR field 
equations for a curved spacetime in the 3 + 1 split under 
the assumption of the conformal flatness condition (CFC) 
for the three-metric. The hydrodynamics equations are 
formulated in conservation form, and are evolved with 
high-resolution shock-capturing schemes. 

In the following subsections, we summarise the mathe- 
matical formulation of the metric and hydrodynamic equa- 
tions, and the numerical methods used for solving them. 



2.1 Metric equations 

We adopt the ADM 3+1 formalism of 
lArnowitt. Deser fc Misnerl (| 19621 ) to foliate a spacetime en- 
dowed with a four-metric into a set of non-intersecting 
spacelike hyp ersurf aces. The line element is then given by 

ds 2 = <? M „ dx^dx v = —a 2 dt 2 +^/ij(dx'+fl'dt)(dx : '+[3 : 'dt), (1) 

where a is the lapse function, /3 l is the spacelike shift three- 
vector, and 7ij is the spatial three-metric. 

In the 3 + 1 formalism, the Einstein equations are split 
into evolution equations for the three-metric 7^ and the ex- 
trinsic curvature Kij, and constraint equations (the Hamil- 
tonian and momentum constraints) which must be fulfilled 
at every spacelike hypersurface. 

The fluid is generally specified by means of the rest- 
mass density p, the four- velocity , and the pressure P, 
with the specific enthalpy defined as h — 1 + e + P/p, where 
e is the specific internal energy. The three-velocity of the 
fluid as measured by an Eulerian observer is given by v 1 = 
u l /(au°) + f3' /a, and the Lorentz factor W = au° satisfies 
1/Vl - viv* 



the relation W 

Based on the ideas o f IsenbereJ (1 19781 ) and 
I Wilson. Mathews fc Marronettil dl996fh an d as d one in 
the work of IDimmelmeier. Font fc Mullerl (|2002al lbh. we 
approximate the general metric <? M „ by replacing the spatial 
three-metric 7^ with the conformally flat three-metric 



(2) 



where jij is the flat metric and 4> is & conformal factor. In 
this CFC approximation, the ADM equations for the space- 
time metric reduce to a set of five coupled elliptic nonlinear 
equations for the metric components, 

A0 = ~2nq> 5 (phW 2 -P)- 5 ^P, 



2 NUMERICAL METHODS 

We construct our initial rotating NS equilibrium models us- 
ing a variant of the self-consistent field method described 



A(q0) = 2na(j) 5 (ph(3W 2 - 2) + 5P) + a<j> 



(3) 



A/3* = IQ-Ka^phW 2 ^ + 20 lo 7r j V,(a 9 r 6 ) - ^\7 i V k f3 k 
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where the maximal slicing condition, K\ = 0, is imposed. 
Here Vj and A are the flat space Nabla and Laplace oper- 
ators, respectively. For the extrinsic curvature we have the 
expression 



0. 



(4) 



which closes the system Q- 

We rewrite the above metric equations in a mathemati- 
cally equivalent form by introducing an auxiliary vector field 
W 1 and obtain 



A<f> = -27T0 5 (phW 2 -P)-<t> 



■K;,l< : ' 



A(a(j>) = 2nacf> 5 (ph(3W 2 - 2) + 5P) + a<j> 



. 7 7K tJ K'- 
8 



(5) 



A/3 1 = 2%{2a<f)- & k i ^ - ^V l V fe /3 fc , 



8tt4> 10 phW 2 v l - Iv l v fe w fc , 



where the flat space extrinsic curvature is given by 



(6) 



and relates to the regular extrinsic curvature as Kij — 4> 2 Kij 
and K ij = 4> 10 K ij . The advantages of this new formulation 
of the metric equations will be discussed in a future publi- 
cation. 

Note that the metric equations do not contain explicit 
time derivatives, and thus the metric is calculated by a fully 
constrained approach, at the cost of neglecting some evo- 
lutionary degrees of freedom in the spacetime metric (e.g. 
dynamical GW degrees of freedom). 

The accuracy of the CFC approximation for isolated 
compact stellar objects has been tested in various works, 
both in the context of stellar core collapse and for equilib- 
rium models of NSs (for a detailed compari son between the 
CFC approximation of GR and full GR. seelOtt et al.ll2007l. 
and r eferences therein). In particular, iDimmelmeier et alj 
(|2006l ) compared collapse simulations of rotating massive 
stellar cores in the CFC approximation w ith the full GR 
calculations of IShibata fc Sekiguchil (|2005l ) . Although such 
massive stellar cores are rather unmotivated astrophysically, 
they are nice toy models where the spacetime dynamics dur- 
ing collapse is violent, similar to what is expected in the case 
of the collapse of NSs to HQSs. For example, some models al- 
most collapse to black holes, with the val ue of the lapse func- 
tion r eaching 0.29. The comparsion by IDimmelmeier et alj 
|2006l ) reveals very good agreement between the CFC and 
full GR calculations. We thus conclude that the spacetime 
of rapidly rotating NS models (whether uniformly or differ- 
entially rotating) is still very well approximated by the CFC 
metric ((2)). The accuracy of the approximation is expected 
to degrade only in extreme cases, such as a rapidly rotating 
black hole, a self-gravitating thin disc or a compact binary 
system. 



0. 



(7) 



where J M = pit M is the rest-mass current, and V,j denotes 
the covari ant derivative with re spect to the four-metric <; M „. 
Following iBanvuls et al.l l|l997t ) , we introduce a set of con- 
served variables in terms of the primitive (physical) variables 
(fi,Vi,e): 



D = pW, S l =phW 2 v 1 , 



r = phW 2 - P-D. 



Using these, the local conservation laws (O can be written 
as a first-order, flux-conservative hyperbolic system of equa- 
tions, 



dt dx 

with the state vector, flux vector, and source vector being 
U = \D,Si,rl 



(8) 



F l = \pv\ SjV + 5}P, tv + Pv 
S = 



o,|t 



dxi 



(9) 



dt 

dxi 



+ T lJ K l 



respectively. Here v l = v l — (5 % /a, and y/—g = a^, with 
g — det(g^ v ) and 7 = det(7ij); the r^ v are the Christoffel 
symbols associated with 

The system of hydrodynamic equations (|HJ is closed by 
an EoS, which relates the pressure to some thermodynami- 
cally independent quantities; in our case P = P(p,e). 



2.3 Numerical methods for solving the metric and 
hydrodynamics equations 

The hydrodynamic solver performs the numerical time inte- 
gration of the system of conservation equations (JSJ using a 
high-resolution shock-capturing (HRSC) scheme on a finite- 
difference grid. In (upwind) HRSC methods, a Riemann 
problem has to be solved at each cell interface, which re- 
quires the reconstruction of the primitive variables (p, v l , e) 
at these interfaces. We use the PPM method for making the 
reconstruction, which yields third-order accuracy in space 
for smooth flows and away from extrema. The numerical 
fluxes are co mputed by means o f Marquina's approximate 
flux formula (|Donat et al.l ["L998T ). The time update of the 
state vector U is done using the method of lines in com- 
bination with a Runge-Kutta scheme having second-order 
accuracy in time. Once the state vector has been updated in 
time, the primitive variables are recovered using an iterative 
Newton-Raphson method. To numerically solve the elliptic 
CFC metric equations (J3J) we utilise an iterative nonlinear 
solver based on spectral methods. The combination of HRSC 
methods for the hydrodynamics and spectral methods for 
the metric equations within a multidi mensional numerical 
code has been presented in detail in IDimmelmeier et al.l 
(|2005l V 



2.2 General relativistic hydrodynamics 

The hydrodynamic evolution of a standard relativistic per- 
fect fluid is determined by the local conservation equations 



1 Note that here we use an analytically equivalent reformulation 
of t he energy source term as compared w ith the one presented in 
e.g. IDimmelmeier. Font fc Mullej ll2002ah . 
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The CoCoNuT code uses Eulerian spherical polar coor- 
dinates {r, 8}; for the models discussed in this work we are 
assuming axisymmetry with respect to the rotation axis and 
also equatorial symmetry. The finite-difference grid consists 
of 200 radial grid points and 40 angular grid points, which 
are equidistantly spaced. A small part of the grid, which 
initially corresponds to 60 radial gridpoints, covers an arti- 
ficial low-density atmosphere, extending beyond the stellar 
surface, whose rest-mass density is 10~ 17 of the initial cen- 
tral rest-mass density of the NS. 

Since the calculation of the spacetime metric is compu- 
tationally expensive, it is updated only once every 25 hydro- 
dynamic timesteps during the evolution and is extrapolated 
in between. The suitability of this procedure has been tested 
and di scussed in detail in iDimmelmeier. Font fc Mullerl 
(2002a). We also note that tests with different grid reso- 
lutions were performed to check that the regular grid reso- 
lution specified above is appropriate for our simulations. A 
check on the relative violation in the conservation of total 
rest mass, ADM mass and total angular momentum showed 
that each of these quantities is conserved to within 1% dur- 
ing the entire evolution time. 



2.4 Gravitational waves 

The GWs emitted by the collapsing NS are com- 
puted using the Newtonian quadrupolc formula in its 
first time-integrated form (the first-moment of mo- 
mentum density formulation as de scribed in detail in 
Dimmelmeier. Font fc Mullerl l2002bT ) in the variant of 
Shibata fc Sekiguchil ( 2004 ). This yields the quadrupole 



wave amplitude A20 as the lowest-order term in a multi- 
pole expan sion of the ra diation field into pure-spin tensor 
harmonics (|Thornelll98"(jh . The wave amplitude is related to 
the dimensionless GW strain h in th e equatorial plane by 
j|Dimmelmeier, Font fc Mullerj [2002bT ) 



h 



iE2 



HAM = g 8524 x 10 -21 

8 v x r \W a cm 



10 kpc 



(10) 



where r is the distance from the emitting source. 

We point out that although the quadrupole formula is 
not gauge invariant and is only strictly valid in the Newto- 
nian slow- motion limit, for GWs emitted by pulsations of ro- 
tating NSs it gives results that agree very well in phase and 
to about 10% - 20% in amplitude with more sophistic ated 
methods (|Shibata fc Sekiguchilliool iNagar et alj|2007l ). 



3 STELLAR MODELS AND TREATMENT OF 
THE PHASE TRANSITION 

3.1 Initial neutron star model 

Following LCCS, we compute the initial equilibrium NS 
model before the phase transition using a polytropic EoS, 



P = Kp< 



(11) 



where K and 7 are constants. We choose 7 = 2 and K — 100 
(in units where Mq = 1) for all of the GR initial models 
considered in the present study. On the initial timeslice, we 
also need to specify the specific internal energy e. For a 



polytropic EoS, the thermodynamically consistent e is given 
by 

K 



„7-l 



7-1 



(12) 



3.2 Hybrid quark star model 

Due to the complexity of the fundamental theory of strong 
interactions, all theoretical studies of quark matter in com- 
pact stars are based on phenomenological models. The MIT 
bag model EoS, whic h has been u s ed extensively for this 
(for a review see, e.g. IWeber 1999]; iGlendennind I2OO2T ) ■ is 
based on the following assumptions: (i) quarks appear in 
colour neutral clusters confined to a finite region of space, 
the volume of which is limited by the negative pressure of 
the QCD vacuum; (ii) within this region, interactions be- 
tween the quarks are weak and can be treated using low- 
order perturbation theory in the coupling constant. These 
two assumptions allow the two main features of QCD to be 
modelled, namely colour confinement and asymptotic free- 
dom. The parameters of the bag model are the bag constant 
B, the masses of the quarks and the running coupling con- 
stant a s , whose value at the scale rel evant for typical qu ark 
chemical potentials is a 3 £ [0.4,0.6] (|Benhar et al.ll2007l ). 

At the moment there is no general consensus about 
the value of B. Fits to the spectrum of li ght hadrons give 
B 1/4 w 145 MeV (|DeGrand et all 1 19751 ). while the ad- 
justment of B wi th hadronic structure functions sug gests 
B 1/A ~ 170 MeV (|Steffens. Holtmann fc Thomas|[l99a ). On 
the other hand, lattic e QCD calc ulations predict values up 
to B 1/4 ~ 190 MeV (|Satj|l982i ). Hereafter, and following 
LCCS, we take B 1/4 = 170 MeV. 

The masses of the u and d quarks are of the order of a 
few MeV (|Hagiwara et al.ll2002l ) and can therefore be mostly 
neglected, whereas the mass of the s quark is much larger, 
its value being in the range m s £ [80, 155] MeV. Neverthe- 
less, including this mass for the s quark, rather than taking 
it to be massless, would decrea se the pressure by only a few 
percent (|Alcock fc O linto 1988]). We do not expect that this 
would change our results qualitatively, and so we neglect it 
in our study. We also neglect the residual interaction be- 
tween the quarks and approximate their temperature as be- 
ing zero. The EoS of the MIT bag model for massless and 
non-interacting quarks at zero temperature is given by 



° q = -(e-4B), 



(13) 



where e is the total energy density. 

A fundamental problem appears to arise in using the 
hydrodynamic equations of Section [2] to describe the quark 
medium since the quarks are being treated as having zero 
rest mass. The quantity e represents only the internal en- 
ergy density of the quarks and contains no rest-mass con- 
tribution. Rest mass appears as a fluid quantity, within this 
picture, only when the quarks become confined. Neverthe- 
less, the continuity equation [the first equation of the system 
of the hydrodynamic equations (JSJ] remains well-defined if 
thought of in terms of baryon number (which is defined 
for the quark medium) rather than in terms of rest mass. 
In order to have a unified treatment for the quark matter 
and the hadronic matter (which is necessary since we have 
transformation between the two), one can define a quantity 
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p = nb?nb in the quark medium (where rib is the baryon 
number density and nib is the rest mass per baryon in the 
hadronic medium) and then formally split the internal en- 
ergy density of the quark medium into two parts, writing 
e = p + pe as usual but bearing in mind that the first term 
on the right hand side represents just a part of the inter- 
nal energy density in the quark phase. If one does this, it is 
easy to show that the treatment of the hydrodynamics goes 
through unchanged in a consistent way, using this p and e. 

For the normal hadron matter, during the evolution we 
use an ideal-gas type of EoS 

Ph = (7 - l)pe. (14) 

However, in contrast to LCCS, we do not reduce the adia- 
batic index 7 from its initial value of 2, because we see no 
physical mechanism which could be responsible for a global 
reduction of pressure in the hadronic matter phase. Conse- 
quently, in our GR models the collapse is caused solely by 
the pressure change due to the introduction of quark matter 

in the core of the NS. 

As first shown bv lGlendenrdnel (|l99ll . [l992T ), if the sur- 
face tension between the phases is not too large, relaxing 
the condition of local electrical neutrality would allow for 
the possibility of having coexistence between the two phases 
within a certain range of densities. In a region where this 
applied, one would then have many inter-mixed microscopic 
zones of the lower-density hadronic matter and of the higher- 
density quark matter. Each zone would have a net electric 
charge but with the mixture being electrically neutral on av- 
erage and with the volume fraction occupied by the higher- 
density phase growing from at the lower-density boundary 
of the mixed phase ph m up to 1 at the upper-density bound- 
ary p qm . 

The value of the lower threshold density ph m , above 
which free quarks start to appear, is rather uncertain. From 
simple geometrical considerations, nucleons should begin to 
touch at p ~ (47rr„ uc /3), which for the characteristic nu- 
cleon radius r nuc ~ 1 fm gives a few times nuclear sat- 
urati on density, p nuc = 2.7 x 10 14 g cm~ 3 (|Glendennind 
1 19891 ). For densities above this, one expects that the bound- 
aries of particles like p, n, E~, A and K~ would dissolve 
and that qua rks would sta rt to populate free states outside 
the hadrons (|Weberill999l ). The value of the upper thresh- 
old density p qm , marking the boundary between the mixed 
phase and the pure deconfined quark matter phase, is also 
uncertain (and is model dependent) but it is very probably 
in the range 4-10p nuc . 

According to this picture, a hybrid star would then be 
composed of either two or three parts: (i) a pure hadronic 
matter phase for p < phm, (h) a mixed phase of the confined- 
deconfined matter for phm < p < p qm , and (iii) a pure quark 
matter phase for p > p qm (this might or might not be present 
in practice, depending on the maximum density reached). 

In this paper, we follow LCCS in adopting the picture 
outlined above and we also follow their prescription in for- 
mulating the EoS for the HQS matter: 

{-Ph for p < p hm , 

oP q + (1 - a)P h for p hm < p < p qm , (15) 

P q for p qm < p, 

where 



Table 1. Summary of the set of initial models: p is the rotation 
period of the NS, Mq is the total rest mass, M is the gravitational 
mass, T is the rotational mass energy, W is the gravitational bind- 
ing energy, p c ; is the central rest-mass density and r e /r p is the 
ratio of the equatorial and polar radii. Note that the initial mod- 
els A5, B4 and C3 are identical. Model N is one of the Newtonian 
initial models used by LCCS. 

Model p M M p c j r e /r p T/\W\ 

[ms] [Af Q ] [M ] [10 14 gcm- 3 ] [%]_ 

Al 1.00 1.98 1.81 11.25 0.635 8.44 

A2 1.20 1.85 1.70 11.25 0.785 5.30 

A3 1.40 1.80 1.65 11.25 0.847 3.75 

A4 1.60 1.77 1.62 11.25 0.885 2.80 

A5 1.80 1.75 1.60 11.25 0.910 2.18 

A6 2.00 1.73 1.59 11.25 0.928 1.74 

A7 2.99 1.70 1.57 11.25 0.968 0.76 

A8 5.98 1.69 1.55 11.25 0.992 0.50 

Bl 1.30 1.75 1.62 8.42 0.746 6.32 

B2 1.40 1.75 1.61 9.48 0.815 4.57 

B3 1.60 1.75 1.61 10.63 0.878 2.98 

B4 1.80 1.75 1.60 11.25 0.910 2.18 

B5 2.00 1.75 1.60 11.80 0.931 1.66 

B6 2.98 1.75 1.60 12.86 0.972 0.66 

CI 1.80 1.65 1.53 8.74 0.882 2.89 

C2 1.80 1.70 1.57 9.88 0.897 2.51 

C3 1.80 1.75 1.60 11.25 0.910 2.18 

C4 1.80 1.80 1.65 13.70 0.927 1.75 

N 1.20 2.20 1.89 9.34 0.695 7.71 



a = 1 - ( Pqm 9 V (16) 

\ P q m — Phm / 

is a factor quantifying the contribution of each of the com- 
ponents of the mixed phase to the total pressure. As stated 
above, we take P q to be given by the MIT bag model JT3J), 
while Ph is calculated using the ideal-gas EoS (|14|l . The pa- 
rameter S is introduced in order to control the quark matter 
contribution to the pressure in the mixed phase: with larger 
S, the contribution from P q increases. For 5 = 1 we recover 
the EoS of LCCS. We again emphasize that, in contrast to 
LCCS, we do not reduce the effective adiabatic index of the 
nuclear matter in our GR models, but rather keep it at its 
initial value 7 = 2 during the evolution. 

For our GR models we define the transition density ph m 
from the pure hadronic matter phase to the mixed phase 
to be where P q vanishes, similarly to LCCS (although they 
were identifying energy density and rest-mass density within 
their Newtonian regime). This corresponds to ph m = 6.97 x 
10 14 g cm -3 = 2.58 pnuc for B 1/4 = 170 MeV and p nuc = 
2.7 x 10 14 gem -3 . Following LCCS, we set p qm = 9 p nU c 
in our simulations (corresponding to 24.3 x 10 14 g cm -3 for 
our value of p nU c), but this is just a rough estimate to give a 
working value. When we make direct comparisons with the 
Newtonian models of LCCS in Section [4] we use their values 
Phm = 7.24 x 10 14 g cm -3 and p qm = 25.2 x 10 14 g cm -3 , 
which slightly differ from our standard values. 

We want to stress that this treatment of the EoS via 
Eqs. (j 151 116|l is extremely rough as a representation of mat- 
ter in the mixed phase experiencing the phase transition, 
especially when one bears in mind the behaviour of fluid el- 
ements undergoing successive compression and decompres- 
sion and changing the proportions of the phases. In partic- 
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ular, it neglects any possible effect from local heating and 
from the creation and subsequent emission of neutrinos dur- 
ing the collapse and the subsequent bounces. However, it 
does have the advantage of being simple and parametris- 
able. By changing the values of the free parameters (e.g. 8, 
Phm or pqm) it allows us to modify easily the properties of 
the EoS. 



3.3 Parameter space 

The properties of our models for the phase-transition- 
induced collapse depend on a number of free parameters 
including the initial rotation period p, the total stellar rest 
mass Mo, the pressure contribution of the quark component 
in the mixed phase (controlled by 8), etc. In order to study 
how the collapse dynamics depends on these quantities, we 
performed simulations for various sequences of models where 
one parameter was held fixed. For instance, in order to inves- 
tigate the impact of rotation we used the models Al to A8 
and Bl to B6 (see Table [TJ. The models of sequence A have 
a fixed central rest-mass density p c ,i = 11.25 x 10 15 g cm~ 3 
and a varying rotation period in the range from p = 1.00 
to 6.18 ms. The ones of sequence B have a fixed rest mass 
Mo = 1.75 Mg and rotation periods from p = 1.30 to 
2.89 ms. The consequences of a variation in the rest mass 
of the initial NS are explored in the sequence of models CI 
to C4, where the rotation period is held fixed at p — 1.80 ms 
and the rest mass takes values from Mo = 1.65 to 1.80 M©. 

These initial models were then evolved as sequences CA, 
CB, and CC with fixed EoS parameters, choosing 8 = 2 in 
each case. For the CD sequence, however, we used different 
values for 8 (from 1 to 5) in order to assess how varying the 
quark contribution to the EoS in the mixed phase influences 
the dynamics. In this sequence, the models CD1 to CD5 use 
the initial model A2, while CD6 to CD9 are based on the 
initial model A5. 

In addition, in order to validate our code and to dis- 
cuss the results obtained by LCCS in more detail, we simu- 
lated some of their models with a Newtonian version of the 
CoCoNuT code. Our Newtonian models CN1, CN2 and CN3 
(all based on the equilibrium model N) are identical to their 
models G1.95, R and G1.75, with the adiabatic index 7 of 
the EoS being reduced everywhere as they had done. For 
these models we used exactly the same EoS as in LCCS, 
which differs from our regular EoS in various ways (see Sec- 
tion I3.2|l . The models labelled CN n i pr and CNi pr are New- 
tonian test models whose properties are discussed in detail 
in Section |4~21 

In Table [2] we list important quantities obtained in the 
simulations of the collapse models which are discussed in the 
following sections. Note that the total evolution time for all 
models is U = 50 ms. At that point we simply terminate the 
evolution, which however could be extended to much longer 
times given the long-term stability of our code. 



c 

o 




4 COMPARISON WITH NEWTONIAN 
MODELS 

Before discussing the results for our GR models, we present 
here our simulations of three of the Newtonian models also 
studied by LCCS. We begin, however, with a description of 



Figure 1. Time evolution of the central rest-mass density p c for 
the Newtonian collapse model CN2 (which is identical to model 
R of LCCS). 



the qualitative features of phase-transition-induced collapse 
of a rotating NS to a HQS which is relevant also for the later 
GR models. 



4.1 Collapse dynamics and gravitational radiation 
waveform 

The EoS of the deconfmed quark matter in the stellar core 
generally gives a smaller pressure contribution than that of 
the hadronic Eo£0 and so the phase transition in the NS 
core leads to an instability of the progenitor NS (which is in 
equilibrium before the transition) and the entire NS starts to 
contract. Depending on the parameters used and the phase 
transition timescale, the infall phase typically lasts between 
0.3 and 0.5 ms. As the pressure in the core rises with in- 
creasing density, the infall decelerates and the contraction 
of the inner core is stopped, while the outer regions are still 
falling in. In the case of a rotating NS, the deceleration of the 
core can be augmented by the increase of centrifugal forces 
due to angular momentum conservation in the contraction 
phase. 

Because of its inertia, the core overshoots its new equi- 
librium configuration, rebounds, expands again and then 
re-collapses. It typically experiences many such distinct se- 
quences of infall, bounce, and re-expansion in the form of 
pronounced, mainly radial ring-down pulsations until it fi- 
nally reaches a new equilibrium state. Fig. [T] shows the time 
evolution of the central rest-mass density p c for the New- 
tonian model CN2, where this oscillatory behaviour can be 
clearly seen. 

If the initial NS is non-rotating, then the newly-formed 
HQS pulsates only radially, and only I = modes are present 
in the oscillation spectrum (unless significant convection de- 
velops). In this case, no GWs are emitted. However, if a 
rotating initial model undergoes a collapse and ring-down, 



2 In principle, the quark matter could give as large a pressure 
contribution as the hadronic matter (see, e.g. lAlford et alj|2007l 
for a recent discussion). However, in our study we do not consider 
those cases. 
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Table 2. Summary of the set of collapse models: 7 is the adiabatic index of the hadronic matter during the evolution, 8 is the EoS 
parameter that specifies the contribution of the quark matter pressure in the mixed phase, p Cj b is the value of the central rest-mass 
density at bounce, |/i| m ax is the maximum value of the GW signal strain during the evolution, E gw is the energy emitted in GWs (during 
a total emission time of tf = 50 ms), fp and faf are the frequencies of the fundamental quasi-radial and quadrupolar modes, respectively, 
while Tp and T2f arc the damping times of those modes in the GW signal. In addition, we give the phases cj>\ and <j>2 of the F-mode and 
the 2 /-mode and the relative amplitude A\/A2 as obtained from a fit to the GW signal according to Eq. ||B1[| . Note that the collapse 
models CA5, CB4 and CC3 are identical, as are models CB2 and CD2, and models CB5 and CD7. During the contraction, models CB6, 
CC4 and CD9 form an apparent horizon and become black holes. The phase-transition-induced collapse models CN1, CN2 and CN3, 
as well as the test collapse models CN n i pr and CNi pr are computed with a Newtonian treatment. Where no values are given for the 
damping times, this signifies that the model either collapses to form a black hole or that no unambiguous damping could be diagnosed 
in the GW signal (mostly due to mode resonance; see Section |5.5|I . 
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then GWs of considerable amplitude can be emitted, as 
shown in the waveform plot in Fig. [2] again for model CN2. 

Comparing a simple collapse model with a purely ideal- 
gas EoS to a regular collapse model with a quark contri- 
bution to the EoS, LCCS demonstrated that the phase- 
transition-induced collapse of a rotating NS to a HQS pre- 
dominantly excites two quasi-normal pulsation modes: the 
fundamental 1 = quasi-radial i^-mode and the fundamen- 
tal I = 2 quadrupolar 2 /-mode. These stand out in the power 
spectrum of many fluid or metric quantities and, in partic- 
ular, in that of the gravitational radiation waveform as pre- 
sented in Fig. [3] All of our collapse models exhibit the pre- 
dominance of these two modes as a generic featur^E For the 

3 We perform the mode identification by perturbing equilibrium 



more slowly rotating models, the contribution from higher 
order modes at higher frequencies can become comparable 
to that of the two fundamental I — 0, 2 modes, but these 
higher frequency modes have damping times that are sig- 
nificantly shorter than those of the fundamental modes and 
their impact on the waveform (and other quantities) dies 
out quickly. 

Comparing our values of selected quantities describing 
the collapse dynamics with the corresponding ones in LCCS 
shows that our code is able to accurately reproduce the orig- 



modcls that have similar structure to those of the collapse rem- 
nants using specific / = and I = 2 trial eigenfunctions and 
analysing the response to these perturbations. This method is 
described in detail in LCCS. 
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Table 3. Comparison of various quantities as calculated in our simulations for the Newtonian collapse models CN1, CN2 and CN3 (top 
row for each model) with the results published by LCCS (bottom row). For each quantity we also give the relative difference between our 
results and theirs. Note that we multiply their values for (h 2 ) 1 / 2 by \f2 to undo the angular averaging and obtain |/i| max . The values of 
p c b for models CN2 and CN3 are extracted from Figs. 6 and 10 in their article, respectively, while for model CN1 they present no data 
from which to read off the central rest-mass density at the bounce. 
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Figure 2. Time evolution of the GW strain h at a distance of 
10 Mpc for the Newtonian collapse model CN2. 
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Figure 3. Power spectrum h (in arbitrary units) of the GW strain 
h for the Newtonian collapse model CN2. The narrow peaks of 
the .F-mode and 2 /-mode clearly dominate the spectrum. 



inal results (see Table [3}, despite being based on a different 
formulation of the hydrodynamic equations and utilising a 
different coordinate system. Furthermore, Figs. [T] [2] and [3] 
which correspond to the data plotted in Figs. 6, 7 (centre 
panel) and 12 (dashed line) of LCCS, also exemplify the 
excellent agreement both qualitatively and quantitatively. 
Making this comparison and demonstrating the good agree- 



ment obtained at a Newtonian level is important for remov- 
ing possible doubts about the analysis which we present in 
the next subsection concerning the onset and development 
of convective instabilities. 



4.2 The role of convection in generating 
differential rotation 

A conspicuous difference from the results of LCCS that we 
observe in the simulations of Newtonian models performed 
with our code is the significantly smaller damping of the 
post-bounce oscillations. This is not only apparent from the 
much larger values that we find for the damping times tf 
and T2f (see Table [3} but can also be noticed by comparing 
our Figs. \T\ and [5] to the corresponding Figs. 6 and 7 (cen- 
tre panel) in LCCS. This has important consequences for 
the physical interpretation given in LCCS for the damping 
of the quasi-radial post-bounce pulsations: they suggested 
that the dominant damping mechanism is conversion of the 
kinetic energy of these pulsations into differential rotation. 
According to their discussion, another less significant part of 
that kinetic energy is lost when matter at the boundary of 
the HQS is ejected by shock waves, while yet another small 
amount of damping is due to numerical dissipation. 

The 3D Newtonian code of LCCS used a coarser grid 
spacing than in our 2D GR calculations and this (together 
with some other possible numerical effects) would have led 
to a higher level of numerical dissipation. If the damping of 
post-bounce pulsations which they saw was indeed mostly 
caused by physical processes such as the transformation of 
kinetic energy into differential rotation or mass shedding at 
the boundary, then we would not have seen much smaller 
damping in our simulations for the same models. We there- 
fore conclude that numerical dissipation actually did play a 
major role for the damping seen in the simulations of LCCS 
(although physical processes certainly also played a role). 

Furthermore, if the observed exponential damping, 
which removes energy from the pulsations at a constant rel- 
ative rate, were directly feeding the generation of differential 
rotation, there should be a simple correlation between the 
evolution time t and the increase in the quantity used as a 
measure of differential rotation by LCCS, 



where f2 



pr sin 8(n- Q) 2 dV, 



(17) 



sin 6 v$ is the local angular velocity with 



10 E. B. Abdikamalov, H. Dimmelmeier, L. Rezzolla and J. C. Miller 



60 



^ 0.5 



o 




Figure 4. Time evolution of the differential rotation measures 
T d (top panel) and A d (bottom panel) for the Newtonian collapse 
model CN2. 



v$ being the rotation velocity of the fluid, and f2 is the 
volume averaged angular velocity of the HQ£0. Note that 
Td does not follow the additive property for energies and 
so is not a fully satisfactory measure of the kinetic energy 
associated with differential rotation. We instead prefer to 
use as our measure of differential rotation, a quantity which 
is a volume- averaged and density-weighted measure of the 
relevant r<j> component of the Newtonian shear tensor, 



A d = 









dr r 



(18) 



The two quantities (|17[) and (|18[) are related and they have 
similar time evolutions, as shown for model CN2 in Fig. U 
(and as found for all of the models considered), but we prefer 
to use Ad because of it having a clearer physical meaning. 

Overall, the evolution of the two measures of differen- 
tial rotation has the following general behaviour: after an 
initial peak associated with the initial collapse and the main 
bounce (with its height correlating with the intensity of the 
bounce), the average differential rotation stays roughly the 
same for around a millisecond (corresponding to several dy- 
namical timescales during which a number of post-bounce 
oscillations occur) following which it grows considerably up 
to a maximum value, then decreases a little and finally os- 
cillates around an almost constant state. This behaviour is 

4 Note that in contrast to LCCS, we averaged Cl at each time 
level and not only once at the time of bounce, because when we 
did the latter, we obtained a very oscillatory behaviour for T d ; 
we are not clear why LCCS did not observe this. 
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Figure 5. Radial profiles of the specific entropy measure K (top 
panel) and of the r-component of its gradient (bottom panel), 
both calculated in the equatorial plane at the initial time for the 
Newtonian models CN2 (solid line), CN n ; pr (dashed line) and 
CN; pr (dash-dotted line). Both the entropy and the radius are 
scaled to give K c = 1 and r e = 1. 



not in accordance with having continuous conversion of pul- 
sational kinetic energy into differential rotation as the main 
damping mechanism. Rather, it suggests that some other 
process is responsible for creating the observed differential 
rotation. 

The initial peak of Td and Ad can be readily explained 
by the fact that any initially uniform rotation profile will 
become non-uniform during the collapse as a result of the 
non-homology of the collapse. On the other hand, the inter- 
mittent behaviour after the initial peak and the saturation 
at a constant value can be interpreted straightforwardly in 
terms of differential rotation caused by large-scale convec- 
tion developing in the HQS several dynamical timescales af- 
ter the initial collapse when it is still pulsating but is already 
close to a new quasi-equilibrium state. 

In rotating stellar models significant convection can oc- 
cur if the Solberg-H0iland sta bility criteria are violated (see, 
e.g. ICerda-Duran et alj 120071. and references therein). For 
our simple EoS, these translate into the condition that con- 
vection can develop if there is a sufficiently strong negative 
radial gradient of specific entropy (depending on the rota- 
tion rate of the HQS). We find that the method of pressure 
reduction employed by LCCS, which involves uniformly low- 
ering 7 in the ideal-gas EoS (|14|l for the hadronic phase 
without adjusting the value for the internal specific energy 
e at the initial time, indeed results in a very large negative 
specific-entropy gradient for their choices of 7. This is shown 
for model CN2 in Fig. [5] (solid line), where we plot the ra- 
dial profile at t = of the (density dependent) measure of 
specific entropy 

p4i - 1) 



K = 



(19) 



in the equatorial plane, assuming an ideal-gas EoS for the 
entire star with a 7 that is reduced from its initial value of 
2 to 1.85. Note that for a polytrope K is identical to the 
polytropic constant K. 

In order to assess unambiguously the occurrence of 
artificially-produced convection and its impact on the devel- 
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Figure 6. Time evolution of the central rest-mass density p c for the Newtonian models CN2 (left panel), CN n j pr (centre panel) and 
CNi pr (right panel). 




Figure 7. Time evolution of the differential rotation measures T<j (top row) and (bottom row) for the Newtonian collapse models 
CN2 (left panels), CN n ; pr (centre panels) and CN; pr (right panels). 



opment of differential rotation, we set up two simple New- 
tonian collapsing test models, which are both based on the 
initial model N and utilise a purely ideal-gas EoS (|14|l , so as 
to simplify the discussion by removing the influence of quark 
matter on the results. In model CN n j pr (with non-isentropic 
pressure reduction) we reduce the pressure by lowering the 
adiabatic index 7 from its initial value to 1.9 (without then 
adjusting e); this creates a strong initial negative specific- 
entropy gradient that is comparable with the one in model 
CN2 (see Fig. [5}. In model CNi pr (for isentropic pressure re- 
duction) 7 remains at its pre-collapse value of 2 during the 



evolution, while the pressure reduction is now realized by a 
lowering the polytropic constant K by 10%, which keeps the 
specific entropy uniform throughout the entire NS. 

Note that the parameters in the EoS for these two test 
models have been chosen in such a way that the collapse 
dynamics (represented by the maximum density reached at 
the first bounce and the amplitude of the post-bounce pulsa- 
tions) is comparable with that of model CN2, as can be seen 
from Fig. [6] Consequently, the gravitational radiation wave- 
forms of these models also have amplitudes and waveforms 
similar to those of model CN2. 
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Figure 8. Radial profile of the rotation velocity in the equatorial plane for the the Newtonian collapse models CN2 (left panel), 
CN n ; pI . (centre panel) and CN; pr (right panel) at the initial time (solid lines) and at t = 10 ms (dashed lines). For models that develop 
considerable convection, the initially uniform rotation profile is locally destroyed. 



However, because of the different behaviour regarding 
convective instability, the dynamics of the three models 
CN2, CN n ip r and CNi pr fall into two very distinct classes, de- 
pending on whether the initial pressure reduction creates a 
strong negative specific-entropy profile or leaves the specific 
entropy constant. In all models, analysis of the meridional 
velocity fields shows no noticeable convection being present 
at early times (e.g. at the time of the main bounce, around 
t — 0.2 ms). In the isentropic collapse model CNi pr con- 
vection continues to remain unimportant also at later times 
and, accordingly, both Td and Ad remain very small at all 
times (see the right-hand panel of Fig. [7J. In stark con- 
trast, the non-isentropic models CN2 and CN n i pr develop 
several convection vortices in the bulk of the collapsed star 
at t ~ 1 ms (corresponding to a few dynamical timescales), 
which is when the two differential rotation measures start 
to increase. This convection grows rapidly and reaches sat- 
uration almost instantaneously, with typical maximum con- 
vection velocities of ~ 0.03 c. During the entire period when 
Td and Ad axe increasing, the convection remains practically 
constant and redistributes angular momentum and entropy 
locally within each vortex. By times t > 5 ms this has led to 
the specific entropy being almost constant within the spatial 
scale of each convection vortex. Convection then subsides 
again, and Td and Ad remain approximately constant from 
then on. It is clear, therefore, that the distinct phases in the 
time evolution of the two quantities reflect very accurately 
the distinct phases of convection, which acts as the mecha- 
nism that redistributes angular momentum and thus creates 
differential rotation. 



The impact of angular momentum redistribution by 
convection on the initially uniform rotation profile can be 
seen from Fig. [H] which shows plots of the rotation veloc- 
ity v,f, in the equatorial plane for the different models. For 
the strongly convective models CN2 (left panel) and CN n i pr 
(centre panel), the curves are driven away from the initial 
uniform rotation (straight-line) profile to reach a step-like 
profile at late times, whereas the essentially non-convective 



model CNi pr (right panel) maintains an approximately uni- 
form rotation profilqj. 

It is worth stressing that the occurrence of convection 
here, and thus the creation of the differential rotation, is es- 
sentially independent of the presence and strength of post- 
bounce pulsations. Indeed, setting aside some small spurious 
convection close to the stellar boundary (caused by interac- 
tion of matter with the artificial low-density atmosphere), 
the strong bulk convection is caused by having a negative 
specific-entropy gradient and can be switched on and off at 
will, depending on whether the initial pressure reduction de- 
stroys isentropy or not. 

In order to demonstrate this connection further, we con- 
structed initial equilibrium models with a local polytropic 
EoS, 

P = K{p) p\ (20) 

with 7 = 2, where the polytropic "constant" K{p) depends 
on the rest-mass density and thus varies inside the NS. By 
adjusting K(p) we can thus create models with negative 
(as well as positive) specific-entropy gradients of arbitrary 
strengths. When these initial models are evolved with an 
ideal-gas EoS (|14|l . they remain essentially in equilibrium, 
pulsating with only very small amplitudes. However, during 
the evolution they develop convection, and subsequently dif- 
ferential rotation, with a strength that directly corresponds 
to the strength of the negative specific-entropy gradient im- 
posed initially. This is illustrated in Fig. [5] where we show 
plots of the meridional velocity field superimposed on the 
magnitude of the ^-component of the vorticity. For this fig- 
ure we selected three such equilibrium models for which the 
initial specific-entropy profile varies between constant spe- 
cific entropy and a gradient that is comparable to model 
CNfl Clearly, convection is practically nonexistent in the 
isentropic model (left panel), whereas the model with mod- 
erate initial non-isentropy (centre panel) develops consider- 

5 This rather good preservation of uniform rotation is a conse- 
quence of the infall being nearly homologous here. 

6 Note that these equilibrium models have central densities com- 
parable with that of CN2 but they are more compact. 
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Figure 9. The velocity field v r g in the meridional plane (with the rotation axis being in the vertical direction) and the colour-coded 
magnitude of the i/j-component of the vorticity are shown for equilibrium models whose initial specific-entropy profile is, respectively, 
constant (left panel), comparable to model CN2 (right panel) and roughly halfway between those two cases (centre panel). The evolution 
time is at t = 3.5 ms when convection (if present) has saturated. The length-scale of the vectors is the same for all three models, with 
the maximum length corresponding to ~ 0.03 c. The vorticity is measured in arbitrary units. 




t [ms] 

Figure 10. Time evolution of the differential rotation measure Ad for equilibrium models whose initial specific-entropy profile is, 
respectively, constant (left panel), comparable to model CN2 (right panel) and roughly halfway between these two cases (centre panel). 
For the isentropic model, Ad remains very close to zero at all times and is barely visible in the plot. 



able convection in those parts of the NS that are not too far 
from the rotation axis. Finally, for the model with the strong 
negative initial specific-entropy gradient (right panel), con- 
spicuous convection vortices occur throughout the NS. Plots 
of the vorticity for collapsing models show similar patterns 
but there the meridional velocity field is also affected by the 
contribution from the large quasi-radial pulsations. 

Also for these equilibrium models the time evolution of 
Ad exhibits the expected behaviour, as shown in Fig. 1101 (we 
here no longer plot Td which, however, exhibits a very similar 
time evolution to that of Ad)- Note that the quasi-periodic 
modulation of Ad is not caused by pulsations of the star, as 
their amplitudes are too small to be visible in Ad and they 
have higher frequencies. Instead, these oscillations (which 
can also be seen in a power spectrum of Ad for the collapse 
models CN2 and CN n i pr ) are due to temporal variations in 
the vortices, with their timescale being determined by the 
typical convection velocity and the average vortex size. 



In summary, we find that the differential rotation re- 
ported by LCCS is almost exclusively caused by the tran- 
sient convection that occurs if a negative specific-entropy 
gradient is generated by the initial pressure reduction. 
Therefore, we are convinced that the conclusion drawn by 
LCCS about the link between the damping of the large am- 
plitude post-bounce pulsations and the creation of differen- 
tial rotation, although seemingly plausible, is not correct. 
They observed that the kinetic energy stored in the pulsa- 
tions is approximately equal to the maximum value of Ad for 
both model CN2 and model CN3 (their models R and G1.75; 
see their Fig. 8) but this is an unfortunate coincidence. Our 
models CN n i pr and CNi pr demonstrate that convection and 
thus the maximum value of Ad can vary enormously for 
roughly constant pulsation amplitude. 
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Figure 11. Radial profile of the initial pressure P\ in the equa- 
torial plane for the GR collapse model CA5 without any pressure 
reduction (solid line), with the regular EoS treatment (combina- 
tion of quark and ideal-gas EoS with 7 = 2; dashed line) and 
with the EoS description used by LCCS (combination of quark 
and ideal-gas EoS with 7 = 1.85; dash-dotted line). In the latter 
case P is reduced throughout the star and not just where p > Phm- 
The vertical dotted lines mark the radius = 4.35 km where 
P = Phm an d the equatorial radius r c = 11.2 km of the initial NS. 

5 GENERAL RELATIVISTIC COLLAPSE 
MODELS 

5.1 Collapse dynamics and gravitational radiation 
waveform 

We next present our results for the GR models, employing a 
quark contribution to the EoS that is slightly different from 
that used for the Newtonian models (see Section f3 . 2 [) . Over- 
all, we have performed simulations for 23 different models, 
and in Table[5]we summarize the most important quantities 
which characterize the dynamics of each model. 

As discussed previously, our prescription for the trig- 
gering of the phase transition and of the subsequent col- 
lapse differs from the one proposed by LCCS in that we do 
not change the adiabatic index 7, leaving the EoS in the 
hadronic phase unchanged. Although still very idealized, we 
believe that this represents a more consistent description 
of the physics of the phase transition. As a consequence of 
this prescription, the phase transition and collapse in our 
GR models is solely caused by the lower pressure exerted by 
matter which is transformed to the quark phase. The differ- 
ence between our approach and that of LCCS is exemplified 
in Fig. 1111 where we show the different prescriptions for the 
initial pressure reduction when applied to the representative 
model CA5. 

With our prescription, only a small (central) part of 
the NS loses pressure support, and thus the NS will not, in 
general, collapse to a black hole even if the initial model is 
close to the stability limit. Rather, it reaches a new stable 
equilibrium state in the form of an HQS. For instance, in 
model CA5 only 0.49 Mq out of a total mass of 1.75 M is 
subject to the pressure reduction. Nevertheless, the change 
of the central rest- mass density during the contraction from 
its initial value p Ci i to p c .b at bounce reaches values of up 
to ~ 50% for some models (see Tables [T] and [2} , which is 
comparable to what was obtained by LCCS in their Newto- 



nian models with a larger overall pressure reduction. There 
are at least two different reasons behind this large and lo- 
cal increase of p: firstly, the stronger gravitational force that 
the NS experiences in a GR framework naturally amplifies 
the strength of the collapse. Secondly, our initial equilib- 
rium models (in particular the ones with high /0 C) i) are al- 
ready close to the limit beyond which the F-mode becomes 
unstable. For these models, therefore, even a moderate per- 
turbation can lead to a strong local contraction and trigger 
post-bounce oscillations of significant amplitude 

Despite this conceptually important difference in the 
way that the phase transition (and hence the resulting col- 
lapse) is triggered, no major qualitative differences appear 
when comparing results from the GR simulations with those 
from the Newtonian simulations of LCCS. This is shown for 
the representative model CA5 in Figs. [T^] and 1131 where we 
plot the time evolution of the central rest-mass density p c 
and the GW strain h, respectively. We point out, however, 
that in our GR models of sequences CA, CB, CC, and CD 
we observe a small spike in the time evolution of the cen- 
tral rest mass density p c at about 0.1 ms (see Fig. [12}. This 
spike is caused by a density compression wave triggered by 
the non-uniform pressure reduction in those models, which 
at the start of the evolution leads to a noticeable gradient 
in the first radial derivative of the pressure at the interface 
between the mixed and pure hadronic matter phases. 

As in the Newtonian case, the waveform is mainly com- 
posed of the fundamental I = quasi-radial F-mo&e and of 
the fundamental 1 = 2 quadrupolar 2 /-mode (see the Fourier 
spectrum of the GW signal in Fig. I14|) . However, in contrast 
to the Newtonian models, the F-mode is now at a lower 
frequency than the 2 /-mode (with the F-mode now having 
frequencies between 0.87 and 1.19 kHz for our selection of 
models, and the 2 /-mode having frequencies between 1.76 
and 2.15 kHz). This difference is a consequence of the dif- 
ferent density profile in the GR case and is in agreement 
with previous i nvestigations of pulsating 7 = 2 equilib rium 
polytropes (see lDimmelmeier. Stergioulas fc Fontll2006l , and 
references therein). The prominent peak next to that for 
the /-mode is a nonlinear self-coupling of the F-mode at 
twice the original frequency, which (like several other such 
nonlinear modes) is strongly excited due to the violent na- 
ture of the collapse. Using the fitting procedure described 
in Appendix Q3j we have extracted from the waveforms the 
damping times for these two modes, obtaining values be- 
tween tf — 8 and 687 ms, and Tit = 18 and 130 ms, respec- 
tively. Because of the much smaller numerical dissipation of 
our code, the damping times computed for the GR models 
are considerably longer than those for the Newtonian models 
calculated by LCCS. 

Another important quantitative difference with respect 
to the Newtonian models appears in the maximum GW 
strain |/i| max that is significantly smaller here for compa- 
rable overall rotation rates. The deeper gravitational poten- 
tial well in GR gives a shorter contraction timescale and 
higher densities at bounce, which tend to amplify the GW 
signal, while this is counteracted by the more compact col- 
lapsed remnant in GR having a smaller quadrupole moment 



7 For a non-rotating HQS with the EoS of Eq. (1 1 5 I t the unstable 
branch starts at p c = 29.6 X 10 14 g cm -3 . 
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Figure 12. Time evolution of the central rest-mass density p c 
for the GR collapse model CA5. 
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Figure 13. Time evolution of the GW strain h at a distance of 
10 Mpc for the GR collapse model CA5. 
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Figure 14. Power spectrum h (in arbitrary units) of the GW 
strain h for the GR collapse model CA5. The narrow peaks of the 
F- and 2 /-modes (along with a prominent nonlinear self-coupling 
of the F-mode) clearly dominate the spectrum. 



(fo r a detailed discussion of these two competing effects, 
see lDimmelmeier. Font fc Mii llcr 2002b) . However, the main 
reason for |/i| m ax here being smaller than for the Newtonian 
models is that within our scenario for the destabilisation, 
a smaller proportion of the total matter content of the NS 
is directly involved in the collapse and in undergoing large 
density variations. For a source at 10 Mpc, |/i| m ax ranges be- 
tween 0.08 and 1.72 x 10 - 23 for all of the models considered 
here, while the energy _E gw emitted in GWs during the first 
50 ms of the evolution ranges between 10~ 6 and 10~ 4 M@c 2 . 
A more detailed discussion about the detectability of these 
sources is presented in Section 

5.2 Influence of rotation, total rest mass and 
composition of the mixed phase 

We next investigate the impact on the collapse dynamics of 
varying the values of the main model parameters: the initial 
rotation period p, the rest mass Mo, and the exponent 5 
used in the mixed phase. 

Rotation. As in any gravitational collapse, the influence 
of rotation is twofold. On the one hand, the rotational 
flattening of the collapsing fluid produces an increase of 
the quadrupole moment, which could potentially lead to 
stronger GW emission. On the other hand, centrifugal forces 
also tend to slow down the collapse and, in cases where they 
are strong enough, they can considerably reduce the time 
variation of the quadrupole moment and hence the GW am- 
plitude. 

For the sequence CA, in which the initial central rest- 
mass density p c ,b is kept constant, the dependence of the 
maximum GW strain |/i| max on the rotation rate is rather 
straightforward to interpret (see the top panel of Fig. ll5l and 
Table [5J. Except for models CA3 and CA4, whose expected 
GW emission is altered by resonance effects (as discussed in 
detail in Section |5~5]) . |/i| m ax increases monotonically with 
increasing rotation (which we here measure in terms of the 
ratio T/|W|, as this quantity turns out to remain almost 
constant throughout the evolution). Since in this sequence 
the values of the central rest-mass density at bounce are all 
close to 16 x 10 14 g cm~ 3 , we conclude that centrifugal forces 
do not significantly retard the collapse for these models and 
so there is no significant associated weakening effect for the 
GW emission. 

The effect of rotation is also investigated in the constant 
rest mass sequence CB, for which we observe a slight initial 
increase of |/i| max with increasing rotation (see the bottom 
panel of Fig. I15[l that is then reversed for very rapid ro- 
tation (model CB5, whose GW emission is again enhanced 
by mode resonance, is an exceptional case). We attribute 
this behaviour to the fact that in this sequence the central 
rest-mass density, both in the initial model and at bounce, 
drops significantly as rotation grows along the sequence, re- 
sulting in a much less violent collapse. This property can be 
seen in the central over-density at bounce, which amounts 
to p c ,b/p c ,i — 1 = 39% for model CB5, but reaches only 13% 
for model CB1. The interplay of growing rotation and de- 
creasing collapse strength then explains the behaviour of the 
GW peak amplitude seen for sequence CB. 

The influence of rotation on the frequencies of the F- 
mode an d the 2 /-mode can be compared d i rectly with the 
results of lDimmelmeier. Stergioulas fc Fontl (|2006T ). who ob- 
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Figure 15. Dependence of |/i| m ax on the initial rotation strength 
T/IWI for the collapse model sequences CA (top panel) and CB 
(bottom panel). Each model is denoted with a filled circle, and 
models with enhanced GW emission due to mode resonance (see 
Section 15.51 1 are shown with an additional large circle. 
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Figure 16. Dependence of the mode frequencies fp, f2j- and 
f2F on the initial rotation strength T/|iy| for the collapse model 
sequences CA (top panel) and CB (bottom panel). Each model is 
marked by a filled circle. 



tained relations between the frequencies of these two modes 
and the strength of rotation for pulsating equilibrium mod- 
els of NSs with a 7 = 2 polytropic EoS, both for sequences 
with constant central rest-mass density and with constant 
rest mass. In agreement with that work, for our dynamical 
collapse models, /f decreases noticeably with more rapid 
rotation for sequence CA and CB, as shown in Fig. \W\ (see 
also Table [2}. Also, for f2f our models reproduce the initial 
increase and subsequent decrease with growing rotation for 
the constant ini tial central rest-mass density sequence CA 
(sequence BU in lDimmelmeier. Stergioulas fc Fontll2006t ) as 
well as the monotonic decrease with increasing rotation 
for the constant rest mass sequence CB (their sequence 
AU). Our results therefore demonstrate that studies of lin- 
ear pulsations of equilibrium models can be used also in 
the more general context of stellar gravitational collapse 
to predict the dependence of the fundamental mode fre- 
quencies on rotation. In our models, the creation of differ- 
ential rotation by non-homologous contraction during the 
collapse phase is rather small. However, even if the devi- 
ation from uniform rotation wer e stro nger, the study of 
iDimmelmeier. Stergioulas fe Font] (|200fj ) suggests that the 
influence on the mode frequencies would still be weak. 

In Fig. [16] we also show the behaviour of the nonlinear 
self-coupling of the F-mode at twice the original frequency 
(the 2 - F-mode) which is strongly excited in all of our mod- 
els (see also Fig. I14[) including the Newtonian ones. Such 
nonlinear couplings of linear quasi-normal modes were also 
observed in the s tudy of initially linearly perturbed equ i- 
librium models bv lDimmelmeier. Stergioulas fc Fontl (|2006f) . 
although at much lower excitation levels. In contrast to the 
initial low-level perturbations used in their models, in our 
case the strong collapse and rebound at bounce manages to 



channel a large amount of energy into this particular mode. 
A detailed discussion of the impact of exciting the 2 -F-mode 
for GW emission is presented in Section \5. 5 1 
Rest mass. For their Newtonian models, LCCS reported that 
along a sequence with constant rotation period p, the maxi- 
mum GW strain |/i| max first increases with Mo and then de- 
creases again. In contrast, we find that for our correspond- 
ing sequence CC, |/i| ma x increases monotonically with Mo 
(see Table [2][f|. This different behaviour is probably due to 
a fundamental difference between our GR models and their 
Newtonian counterparts. Besides the obvious differences in 
the structure of the HQS and hence in the mode frequen- 
cies, one should bear in mind that the rest mass of the GR 
models has an upper limit corresponding to the point of 
transition to the unstable branch of equilibrium solutions. 
This is particularly evident for sequence CC, where model 
CC4 does not reach a stable configuration after the collapse 
but instead produces a black hole as deduced from the ap- 
pearance of an apparent horizon. As a result, the eventual 
decrease of |fe| ma x with increasing Mo in the Newtonian case 
may not occur in GR simply because it might require rest 
masses above the upper limit. 

Composition of the mixed phase. The impact on the collapse 
dynamics of varying the parameter S (and hence varying the 
pressure reduction due to the presence of the quarks in the 
mixed phase) is quite obvious. Both for the collapse models 
CD1 to CD5, which are based on the initial model B2, and 
for the models CD6 to CD9, which use the initial model B5, 
the rest-mass density at bounce p c .b and the amplitude of 



8 We have found the same qualitative behaviour also for another 
sequence with constant period p = 1.4 ms and varying rest mass, 
for which we do not present the data here. 
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the post-bounce oscillations grow when the pressure reduc- 
tion is enlarged by increasing the exponent <5 (see Table [2}. 
Since within each of the two subsequences the initial rota- 
tion rate is constant, the growth of p c ,t>, and hence of the 
compactness, with a greater pressure reduction is directly re- 
flected in stronger GW emission and higher values for |/i| max - 
Also in this case, an exception arises with model CD7 which 
is identical to the already mentioned model CB5. Interest- 
ingly, when 8 = 4, the pressure reduction in the mixed phase 
is so great that the corresponding model CD9 collapses di- 
rectly to a black hole without experiencing a bounce. 

Finally, we note that the F-mode frequencies in the se- 
quence of models CD1 to CD5 first decrease with increasing 
S and then increase again. This may be the result of a near 
cancellation of opposite effects: while a larger pressure re- 
duction produces a higher overall density of the post-collapse 
HQS (potentially resulting in higher values of fp), it also 
leads to more rapid rotation due the comparatively greater 
compactness (which should lower fp ) . We also note that the 
large amplitude quasi-radial oscillations in the last model(s) 
of each subsequence of sequence CD are so strongly damped 
shortly after the collapse, that making a precise determina- 
tion of the frequency fp is difficult. 

5.3 Damping of the stellar pulsations 

Both from the values for the damping time Tp of the fun- 
damental quasi-radial F-mode in the GW signal given in 
Table [2] and from Fig. 1171 where we plot the time evolution 
of p c for the most and least rapidly rotating collapse models 
of sequence CA, it is evident that there can be significant 
damping of the post-bounce pulsations in the HQS, for some 
values of the model parameters. In Table U we also report 
the values for the damping time rp, Pc extracted from the 
time evolution of the central rest-mass density p c : these val- 
ues range from 8 to 200 ms for the models considered here. 
With the exception of those models for which the GW emis- 
sion is dominated by mode resonance, the two estimated 
timescales rp and tf, Pc are very similar. 

In Section 3J we have argued that the kinetic energy 
stored in the quasi-radial post-bounce pulsations is not re- 
sponsible for generating differential rotation, and therefore 
the observed damping cannot be attributed to this mech- 
anism as previously suggested. Since, also, the numerical 
dissipation of our code is much too small to be the main 
cause for the strong damping observed in some models (and 
anyway affects all models in the same way) it is necessary 
to offer an alternative physical explanation for why the pul- 
sations are strongly damped for some models, with rp tPc 
being a few ms (and thus comparable with the dynami- 
cal timescale), while for other models the damping is much 
slower, with tf, Pc ~ 100 ms (and thus orders of magnitude 
longer than both the dynamical timescale and the pulsation 
period of the star; in this case, the damping seems to be lim- 
ited essentially to that caused by numerical viscosity) . Since 
the gravitational radiation back-reaction on the system is 
also negligible (and is not taken into account by our confor- 
mally flat approximation of the metric equations anyway), 
the strong damping seen can only be due to hydrodynamic 
effects. 

A careful analysis of our model parameters reveals that 
damping is significant only for those models that both rotate 
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Figure 18. Location of all GR collapse models in the (p/pk)~ 
(Pc,b/Pc,i — 1) plane. The damping time rp „ c of the F-mode, 
as calculated from the time evolution of the central rest-mass 
density, is encoded in the grey scale, with dark/light grey corre- 
sponding to strong/weak damping. Only models with very rapid 
rotation (close to the Kepler limit) and violent post-bounce os- 
cillations due to a very dynamical collapse (upper left corner of 
the plot) have short damping times. 



very rapidly and experience strong pulsations. This can be 
seen in Fig. 1181 where the value of tf, Pc for each model is 
given by the grey scale. Only models located in the upper 
left part of the plane spanned by the rotation period ratio 
p/pK and the over-density at bounce p c ,b/pc,i — 1, have the 
very short damping times. 

For models rotating close to the mass-shedding limit 
(the Kepler limit), the effective gravity near to the equator 
is significantly weakened, vanishing at the mass-shedding 
point. As one goes to models having larger amplitude quasi- 
radial post-bounce pulsations, an increasing amount of mat- 
ter is ejected from low latitudes on the stellar surface during 
each oscillation. This matter goes into the initially artificial 
very-low-density atmosphere and creates an expanding en- 
velope of weakly-bound (or even unbound) material around 
the HQS. This mass shedding causes strong damping of the 
pulsations, as pulsational kinetic energy is converted into 
gravitational potential energy of the ejected matter. Since 
polar perturbations (mainly radial ones) are coupled to axial 
perturbations in rotating stars, the damping rapidly affects 
all modes. 

This mass-shedding-induced damping of oscillations in 
rapidly rotating NSs nea r to the Kepler limit was first 
observ ed and discussed by IStergioulas. Apostolatos fc Fontl 
l)2004l ) for weakly pulsating equilibrium polytropic NS mod- 
els in uniform rotation, treated within the Cowling approx- 
imation (neglecti ng the dynamics of the spacetime ) . In a 
subsequent study, iDimmelmeier. Stergioulas fc Fontl (|2006t l 
found that if the Cowling approximation is abandoned and 
the spacetime is dynamically evolved, the effect of mass- 
shedding is drastically reduced, giving a consequent decrease 
in the damping of the pulsations. However, both of these 
studies treated small-amplitude pulsations of equilibrium 
models (working within the linear regime), whereas in our 
collapse models the pulsations can have very large ampli- 
tudes. 

Considering the example of model CD8: this has a large 
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t [ms] 

Figure 17. Time evolution of the central rest-mass density p c for models CA1 and CA8, the most and the least rapidly rotating models 
of sequence CA. The dotted line marks the pulsation damping time tf „ c for model CA1. 



initial pulsation amplitude, with relative variations in p c 
of around 55% (which is about half of the over-density at 
bounce) but it is rotating too slowly to be close to the Kepler 
limit (with p/pk — 3) and does not show strong damping. 
The five models with the shortest damping timescales, how- 
ever, (jF,p c < if = 50 ms) all have both quite large initial 
pulsation amplitudes > 18% and are also rapidly rotating, 
with rotation-period ratios p/pn < 1.75. This is consistent 
with equatorial mass-shedding being the predominant mech- 
anism for the strong damping of the post-bounce pulsations 
seen for models that both rotate close to the Kepler limit 
and pulsate with large amplitude. For models that are not 
affected by equatorial mass shedding, other damping mecha- 
nisms predominate (such as the conversion of kinetic energy 
into internal energy by shocks, nonlinear coupling of modes, 
numerical or physical dissipation) but all of these operate 
on timescales much longer than the dynamical timescakfl 



5.4 The role of convection 

As discussed above, the damping of post-bounce pulsations 
and the development of differential rotation seen by LCCS 
seems to have been mainly a manifestation of the convec- 
tive motion artificially produced in their simulations by the 
way in which they induced the collapse. Real physical con- 
vection might well be induced by energy input coming from 
the phase transition but this has no connection with the 
artificially-induced convection. Calculating the real convec- 
tion is beyond the scope of the present simplified treatments 
and remains a topic for future work. While we have elimi- 
nated from our calculations the main source of artificial con- 
vection present in the LCCS simulations, some artificially- 
induced convection still remains. There are at least two dif- 
ferent origins for this which apply for two different classes 
of initial models. 



9 As far as viscosity is concerned, in the quark phase the shear 
viscous damping timescale is comparable to that of normal nu- 
clear matter, which is of the order of 10 s s for a typical NS, 
while the bulk viscous damping timescale is of the order of sec- 
onds, assuming a rotation per iod p = 1 ms and an s-quark mass 
m a = 100 MeV llMadsenll200ol) . 
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Figure 19. Radial profile of the specific-entropy measure K (top 
panel) and its radial derivative (bottom panel) in the equatorial 
plane at the initial time for the GR model CA1 (solid line), and 
the corresponding profiles for the Newtonian models CN1 (dashed 
line) and CN3 (dash-dotted line). The specific-entropy measure 
and the radius are scaled to give K c = 1 and r e = 1. 



The first origin is related to the CFC approximation it- 
self: when the data from the initial-model solver are mapped 
to the evolution code and the initial- value problem is subse- 
quently re-solved to satisfy the constraints Q, small errors 
due to the CFC approximation create spurious departures 
from constant entropy. This side effect of the CFC approx- 
imation is clearly larger for rapidly rotating and very com- 
pact models such as CA1, and becomes very small for slowly 
rotating models, vanishing in the non-rotating limit. Fig. 1191 
shows the initial radial profile of the specific-entropy mea- 
sure K [cf. Eq. (|19[>] for model CA1, together with those for 
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Figure 20. Time evolution of the differential rotation measure A<j for the GR collapse model CA1 (left panel) and the Newtonian 
collapse models CN1 (centre panel) and CN3 (right panel). 



the Newtonian models CN1 and CN3 discussed previously. 
It can be seen that the CFC approximation produces a neg- 
ative specific-entropy gradient for model CA1 at the start of 
the evolution. While this gradient is very small in the bulk 
of the star (for r/r c < 0.8), near to the surface it becomes 
comparable to (or even greater than) that for model CN1, 
for which 7 was initially decreased from 2 to 1.95. However, 
it is still much smaller than that occurring when reducing 7 
to 1.75 as in model CN3. 

The CFC-induced violation of the Solberg-H0iland cri- 
teria drives this model (and other rapidly rotating ones) to 
develop significant convection, which we measure in terms 
of the general relativistic equivalent of the averaged shear 
[cf. Eq. dHJ], 



dr 



r 



dV, 



(21) 



where — V V3V 3 is the rotation velocity of the fluipF"!. 
Fig.[20]shows the time evolution of Ad for model CA1 (whose 
convection is, in fact, among the strongest for the GR mod- 
els; left-hand panel), as well as for the Newtonian models 
CN1 (centre panel) and CN3 (whose convection is among 
the strongest for the Newtonian models; right-hand panel). 
It can be seen that both the development of the convection 
with time for model CA1 and its maximum saturation level 
are comparable with those for the Newtonian model CN3, 
while the convection in model CN1 is stronger and occurs on 
a much shorter timescale, which is natural since the initial 
non-isentropy is larger in that case. 

The second origin of spurious convection appears in 
slowly rotating models; this is not due to the CFC approxi- 
mation since the errors coming from that are essentially neg- 
ligible in these cases. Instead, convection is produced here 
by the small violation of isentropy caused by interaction 
of the matter just inside the surface of the HQS with the 
surrounding artificial atmosphere. This perturbation slowly 
propagates inwards and finally affects the entire star at late 



10 As already noted for the Newtonian models, the Newtonian 
equivalent of Ad has a qualitatively very similar behaviour to 
that of Td but only the first of these has a clearly-defined physical 
meaning. 



times. Since in uniformly rotating configurations, centrifugal 
forces tend to have a stabilising effect within the Solberg- 
H0iland criteria, it is easy to understand why this process 
occurs only for the slowly rotating models. An example of 
this mechanism in action is given by model CA8, which is 
the slowest rotator in our set of GR models: here Ad reaches 
a maximum that is close to the saturation level of models 
CA1 and CN1 but is still much smaller than the value for 
the strongly convective model CN3. 

In models that are neither very rapidly nor very slowly 
rotating (e.g. model CC1), convection is almost absent and 
Ad remains at very low values, which can be explained as 
resulting from small deviations away from homology during 
the initial contraction and the subsequent post-bounce os- 
cillations. We repeat, however, that consistent treatment of 
other possible real sources of convection remains a major 
topic for future work. 



5.5 Enhanced emission of gravitational waves via 
mode resonance 



In Section 15.21 we have discussed the influence of the initial 
rotation speed on the frequencies of the post-bounce oscil- 
lations and also commented that for models with the same 
rest mass, those with higher rotation generally have higher 
GW strain due to their increased quadrupole moments. This 
is indeed what is seen, for instance, in the sequence CB as 
reported in Table [2] This general behaviour, however, has 
a notable exception in the case of the comparatively slowly 
rotating model CB5, which has by far the largest GW strain 
for any of those in the sequence CB, with the energy emit- 
ted in GWs in the first 50 ms being at least an order of 
magnitude larger than for any of the other models in this 
sequence. 

In the top panels of Fig.[5T]we plot the time evolution of 
the central rest-mass density p c for models CB3, CB4 and 
CB5. It can be seen that the collapse dynamics of model 
CB5 are not qualitatively different from those of the more 
rapidly rotating members of the same sequence. As rotation 
decreases from CB3 to CB5, the strength of the contrac- 
tion increases, leading to higher central over densities p Cj b at 
the bounce and stronger post-bounce oscillations. This, how- 
ever, cannot account satisfactorily for the very large maxi- 
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Figure 21. Time evolution of the central rest-mass density p c (top panels) and the GW strain /lata distance of 10 Mpc (bottom panels) 
for the GR collapse models CB3 (left panels), CB4 (centre panels) and CB5 (right panels). The enhanced GW emission for model CN5 
due to mode resonance is clearly visible. 




Figure 22. Logarithmic power spectrum (in arbitrary units) of the central rest-mass density p c (top panels) and the GW strain h 
(bottom panels) for the same models as in Fig. 1211 The resonance between the 2 /-mode and the 2 -F-mode sets in as their frequencies 
(marked by dotted lines) approach each other. 



mum GW strain amplitude |/i| max of model CB5, nor for the 
growth of h during the first 20 ms after the bounce in this 
case (see bottom panels of Fig. I21|l . 

Rather, the time evolution of h in model CB5 (the 
delayed growth, the saturation and the subsequent decay) 
suggests the presence of a resonance effect among sev- 
eral modes, at least one of which should be an efficient 



emitter of GWs. One of the obstacles to having strong 
GW emission is that the modes which are most strongly 
excited during the bounce are (quasi-)radial modes which 
are not efficient GW emitters; indeed they would not emit 
at all if rotation were not present to introduce non-radial 
contribut ions. However , it was suggested many years ago 
(see, e.g. lEardleyl 1 19851 ) that the pulsational energy from 
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Table 4. Relevant quantities for mass-shcdding-induced damp- 
ing: p is the rotation period of the NS, pk is the rotation period 
of a free particle on a circular orbit at the equator (Kepler limit), 
p/pK is the ratio of these two periods, p c b/Pc.i — 1 is the central 
over-density at bounce compared to the initial value, and rp „ 
and Tp are the F-mode damping times extracted from the time 
evolution of p c and from the GW signal, respectively. Note that 
we omit here those models that collapse to a black hole. For com- 
pleteness, we also include the Newtonian collapse models CN1, 
CN2 and CN3. 



Table 5. The frequencies fij and f%.p of the fundamental 
quadrupolar 2 /-mode and the self-coupling of the F-mode, for 
the collapse models of sequence CB. A/ rc i = f2-F / f'-if — 1 is the 
relative difference between the frequencies of these two modes. 







p 


Pk 


p/pk 


Pc,b/Pc,i ~ 1 


T F,p c 


T F 




H 


[ms] 




[/oj 


[ms] 


[ms] 


CA1 


1 


00 





89 


1.12 


11 


32 


40 


CA2 


1 


20 





75 


1.60 


43 


36 


49 


CA3 


1 


40 





72 


1.94 


41 


93 




CA4 


1 
1 


ou 


n 
U 


/ U 


2.29 


42 


155 


319 


PAc: 


1 

1 


ou 


A 

u 


(;o 
\))J 


2 61 


49 


1 fin 


418 


CA6 


2 


00 





68 


2.94 


42 


150 


270 


CA7 


2 


99 





67 


4.46 


43 


152 


711 


CA8 


5 


98 





66 


9.06 


43 


132 




CB1 


1 


30 





89 


1.46 


13 


113 


99 


CB2 


1 


40 





80 


1.75 


23 


128 


133 


CB3 


1 


60 





73 


2.19 


35 


156 


196 


CB4 


1 


80 





69 


2.61 


12 


160 


418 


CB5 


2 


00 





67 


2.99 


52 


103 


687 


CC1 


1 


80 





79 


2.28 


18 


136 


143 


CC2 


1 


80 





71 


2.43 


29 


62 


71 


CC3 


1 


80 





69 


2.61 


12 


160 


418 


CD1 


1 


40 





80 


1.75 


10 


200 


150 


CD2 


1 


40 





80 


1.75 


23 


128 


133 


CD3 


1 


10 





80 


1.75 


36 


49 


54 


CD4 


1 


10 





80 


1.75 


19 


19 


19 


CD5 


1 


40 





80 


1.75 


62 


8 


8 


CD6 


1 


40 





80 


1.75 


22 


191 


248 


CD7 


2 


00 





67 


2.99 


52 


103 


687 


CD8 


2 


00 





67 


2.99 


110 


53 




CN1 


1 


20 





86 


1.40 


11 


91 


73 


CN2 


1 


20 





86 


1.40 


32 


16 


9 


CN3 


1 


20 





86 


1.40 


52 


1 


1 



a quasi-radial mode, which contains a significant amount 
of kinetic energy but radiates GWs only weakly, could be 
transferred to a much more strongly radiating quadrupolar 
mode by means of re sonance effects and even para- 
metric instabilities. iDimmelmeier. Stergioulas fc Fontl 
j|200fj). as well as iPassamonti et al.l and 



IPassamonti. Stergioulas fc Nagarl ( 20071 ) have discussed 



this possibility in the context of nonlinear coupling of 
quasi-normal modes for nearly-equilibrium models of NSs. 

In order to investigate whether this mechanism might 
be responsible for the enhanced GW emission observed for 
model CB5, we performed a mode analysis for models CB3, 
CB4 and CB5. The power spectrum of the time evolution of 
the central rest- mass density p c , shown in the top panels of 
Fig. 1221 indicates that there is a lot of energy in the 2F-mode 
which has a peak about an order of magnitude higher than 
the corresponding one for the 2 /-mode, located at slightly 
lower frequencies. As rotation decreases from model CB3 to 
model CB5, the frequencies of these two modes get closer 
until the two peaks almost merge for model CB5, with the 
relative difference between the two frequencies decreasing 
to about 4% (see the bottom panel of Fig. [16] and also Ta- 
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Under these conditions of resonance, the 2 -F- mode 
is able to transfer a considerable amount of energy into the 
2 /-mode, as can be clearly seen in the power spectrum of the 
waveform amplitude shown in the lower panels of Fig. 1221 
If the spectra for model CB5 from Fig. \W\ are produced for 
several time windows, this energy transfer between the two 
modes becomes visible. While the peaks of both p c and h cor- 
responding to the 2 /-mode first gradually grow in the early 
phases of the evolution (showing the initial amplification at 
the expense of the 2 -F-mode) and start descreasing only at 
later times, the peaks of the 2 - F-mode always descrease. 

The strong dependence of the maximum GW strain 
|/i|max on the rotation rate, as seen in Table [5] as well as in 
the bottom panel of Fig. 1151 indicates that the resonant be- 
haviour becomes important only when the frequencies of the 
2 /-mode and the 2 -F-mode are very close to each other. This 
small difference between the two frequencies is also respon- 
sible for the clear beating pattern seen in the waveform for 
model CB5 (see the bottom right panel in Fig. I21|l . Fig. [15] 
also shows that while the pulsation energy contained in the 
F-mode (represented by the spectrum of p c in the top panel) 
is always larger than the corresponding one in the 2 -F-mode, 
when more resonance between the 2 ■ F-mode and the 2 f- 
mode is at work, the 2 F-mode becomes a very efficient emit- 
ter of GWs as well (see the spectrum of h in the lower panel), 
surpassing the F-mode here. This is most likely a conse- 
quence of the altering of the 2 • F- mode's previously quasi- 
radial eigenfuncti on by the interaction with the quadrupol ar 
2 /-mode (see, e.g. IDimmelmeier. Stergioulas fc Fontll2006l ). 

Applying this analysis to the complete set of investi- 
gated models, we find that the same resonance between the 
2 /-mode and the 2F-mode is also at work in models CA3 and 
CA4 (this is clearly visible in the top panel of Fig. ll5Jl as well 
as in model CD8. Also in these cases, the waveforms show 
an initial growth in amplitude over several cycles and then 
a strong beating at later times. In the case of model CD8, 
the growth in h due to the mode resonance becomes promi- 
nent long after the bounce, at t > 20 ms, not reaching its 
maximum |/i| max until t ~ 48 ms, when the quasi-radial os- 
cillations have already been damped to around 38% of their 
initial amplitude. In contrast, for models whose GW emis- 
sion is not influenced by this resonance, |/i| ma x is reached 
at the time of bounce or very close to it. In addition, here 

11 We have also performed simulations for models with other 
parameter values in the close vicinity of those for model CB5 and 
found that CB5 actually exhibits almost the maximum possible 
resonance between the 2 - F-mode and the 2 /-mode. 
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the spectral power of both the 2 -F-mode and the 2 /-mode 
decreases at all times if successive spectra are obtained with 
the technique of shifted time windows. 

5.6 Detectability prospects for the gravitational 
wave emission 

Although the maximum GW strain |/i| max in the waveforms 
of our GR models is about an order of magnitude smaller 
than that computed by LCCS for Newtonian models, the 
long quasi-periodic GW emission that is possible for phase- 
transition-induced collapse may still make this scenario a 
plausible source for GW detectors. Assuming that the strong 
post-bounce oscillations are not damped by any other phys- 
ical mechanisms apart from dissipation of kinetic energy by 
shocks, the damping times that we obtain suggest that the 
effective total emission time for GWs can be much longer 
than the time for which we have followed the evolution of 
our models (if = 50 ms), extending to hundreds of oscilla- 
tions. 

The main GW emission modes are the F-mode and 
the 2 /-mode, which have comparable energies in the power 
spectra, except where the resonance between the 2 /-mode 
and the 2 • F-mode is important, when the combined 2 f- 
mode/2 • F-mode dominates the GW spectrum. Model CA5, 
for which we show in Fig. [23]the power spectrum of the time 
evolution of p c (top panel) and of the GW signal (bottom 
panel), is an example where the frequencies of 2 /-mode and 
2 ■ F-mode are already close and the two modes start to 
merge. 

Exploiting the fact that, for most models, the waveform 
can be very well approximated as a combination of two es- 
sentially monochromatic damped sinusoids, it is straightfor- 
ward to construct a phenomenological waveform expressed 
in terms of six parameters: the mode frequencies /f and 
f2f, the damping times tf and T2f, and the initial phases 
and relative amplitude of the two modes (all of which can be 
found from Table [2}. This can then be applied as a template 
in matched filtering data-analysis algorithms so as to search 
for the waveforms in the data stream more effectively. 

In order to assess the prospects for detection by cur- 
rent and planned interferometric detectors, we next calcu- 
late character istic quantities for the GW signal following 
iThornd i|l987l ). Making a Fourier transform of the dimen- 
sionless GW strain h, 

/oo 
e 2n ' ft hdt, (22) 
-oo 

we can compute the (detector dependent) integrated char- 
acteristic frequency 

and the dimensionless integrated characteristic strain 

where Sh is the power spectral density of the detector and 
She = Sh(fc). We approximate the average (h 2 ) over ran- 
domly distributed angles by h 2 , assuming optimal orien- 
tation of the detector. From Eqs. (|23l I24[) the signal-to- 



noise ratio (SNR) can be calculated as h c /[h rnls (f c )], where 
ft rms = V fSh is the value of the rms strain noise for the 
detector (which gives the theoretical sensitivity window). 

In Table [6] we summarize the values of / c , h c and the 
SNR for all of the models (except those which collapse to a 
black hole) for the currently operating LIGO and VIRGO 
detectors and for the future advanced LIGO detector. For 
all of the detectors we consider a source inside our own 
Galaxy at a reference distance of 10 kpc. The proportion 
of NSs that undergo phase-transition-induced collapse at 
some stage in their lifetimes is not well-known. The phe- 
nomenon could occur at, or soon after, the formation stage 
(giving an event rate roughly proportional to that for core 
collapse supernovae) or it could come at a later stage when 
an old NS is spun up and has its mass increased by ac- 
cretion from a binary partner in an LMXB (an interesting 
case but with an event rate which is thought to be very 
much lower, prob ably ~ 10 -5 yr _1 for the Milky Way; see 
IPfahl et al.ll2003T h Another point is that a phase-transition- 
induced collapse occurring for a NS which is not rapidly 
rotating would not be so interesting for our purposes. Even 
under the most extremely optimistic assumption that the 
phase-transition-induced collapse rate equals that for core- 
collapse supernovae, that would only give a rate of up to 
1 per 20 yr for our Galaxy, which is prohibitively low. On 
the other hand, if such an event did occur in our Galaxy, 
for current interferometric detectors of the LIGO class and 
assuming an emission time if = 50 ms, all of our models 
except one have a SNR above 1. For the strongest emitting 
model CB5, where mode resonance significantly enhances 
the GW emission, the SNR even exceeds 10. With the ad- 
vanced LIGO detector, the SNR lies comfortably above 10 
for all models and reaches almost 400 in model CB5. 

For substantially increasing the event rate, it would be 
necessary for the detector to be sensitive to signals com- 
ing from distances out to the Virgo cluster, at ~ 20 Mpc, 
(for which the supernova rate would rise to more than 1 
per year) . However, at this distance the SNR for our models 
drops to well below 1 even for advanced LIGO. Therefore, 
as for GW signals from supernova c ore collapse (see the dis- 
cussion in iDimmelmeier et al.|[2005h . the second generation 
detectors will improve the SNR of a local event, but will not 
increase the event rate much on account of the inhomoge- 
neous galaxy distribution in the local region of the Universe. 
Only third generation detectors will have the required sen- 
sitivity in the kHz range to achieve a robust SNR at the 
distance of the Virgo cluster. 

Note that for most of the models the integrated char- 
acteristic frequency / c given in Table [6] is not very close to 
either of the two main GW emission frequencies /f and /2* . 
This is because / c reflects the frequency dependence of the 
sensitivity of the detector, because nonlinear mode couplings 
and higher order linear modes also contribute to the GW 
signal (although with lower relative amplitudes; see Fig. I23|) 
and, most importantly, because for many models the GW 
power spectrum of the signal exhibits nearly equally strong 
peaks in the F-mode and the 2 /-mode. 

The detector dependence of f c and h c is also illustrated 
in Fig. 1241 where the locations of the GW signals for all of 
the models are plotted relative to the rms strain noise /i rms 
of the current LIGO and VIRGO detectors (left and centre 
panels, respectively) and the future advanced LIGO detector 
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Figure 23. Logarithmic power spectrum (in arbitrary units) of the central rest-mass density p c (top panel) and the GW strain h (bottom 
panel) for the GR collapse model CA5. The linear quasi-normal modes and the nonlinear (self-)couplings are marked with dotted lines 
and labelled. Hi denotes the first overtone of the _F-mode, while 2 pi and 4 pi are the first overtones of the fundamental 1 = 2 and I = 4 
modes, respectively. 



Table 6. Detection prospects for the GWs: f c is the characteristic frequency, h c is the integrated characteristic GW signal strain and 
SNR is the signal-to-noise ratio, each given for the current LIGO and VIRGO detectors, and for the future advanced LIGO detector. 
The values given for all of the quantities assume a total emission time of tf = 50 ms and are dependent on the rms strain noise h TulB of 
the detector. Note that models CB6, CC4 and CD9 collapse to form a black hole as a result of the phase transition and so are omitted 
here. 



Model 


/c,LIGO 


/c.VIRGO 


/c,adv. LIGO 


ftc.LIGO 


^c.VIRGO 


h c , adv. LIGO SNRlIGO 


SNRvirgo 


SNR adv . LIGO 




[kHz] 


[kHz] 


[kHz] 


" 1Q -20 




" 10 -20 




10 -20 


[at 10 kpc] 


[at 10 kpc] 


[at 10 kpc] 










at 10 kpc 




at 10 kpc 




at 10 kpc 








CA1 


0.954 


0.978 


1.029 


2.55 


2.59 


2.68 


7.5 


12.0 


179 


CA2 


1.755 


1.798 


1.864 


5.65 


5.75 


5.92 


6.8 


11.8 


201 


CA3 


1.845 


1.876 


1.923 


8.52 


8.63 


8.80 


9.4 


16.6 


287 


CA4 


1.619 


1.660 


1.735 


3.56 


3.62 


3.75 


4.8 


8.3 


139 


CA5 


1.349 


1.382 


1.450 


1.78 


1.81 


1.87 


3.1 


5.3 


86 


CA6 


1.415 


1.452 


1.527 


1.55 


1.58 


1.63 


2.5 


1.1 


71 


CA7 


1.449 


1.485 


1.572 


0.72 


0.73 


0.76 


1.1 


2.0 


32 


CAS 


1.609 


1.672 


1.788 


0.29 


0.30 


0.31 


0.4 


0.7 


11 


CB1 


1.349 


1.374 


1.417 


1.31 


1.33 


1.36 


2.3 


1.0 


64 


CB2 


1.304 


1.335 


1.388 


2.12 


2.16 


2.22 


3.9 


6.7 


108 


CB3 


1.322 


1.353 


1.411 


2.09 


2.12 


2.18 


3.8 


6.1 


101 


CB4 


1.349 


1.382 


1.450 


1.78 


1.81 


1.87 


3.1 


5.3 


86 


CB5 


2.073 


2.078 


2.086 


13.18 


13.27 


13.24 


12.3 


21.9 


389 


CC1 


1.370 


1.394 


1.435 


1.06 


1.08 


1.10 


1.8 


3.1 


51 


CC2 


1.370 


1.394 


1.442 


1.48 


1.50 


1.54 


2.6 


1.1 


71 


CC3 


1.349 


1.382 


1.450 


1.78 


1.81 


1.87 


3.1 


5.3 


86 


CD1 


1.346 


1.370 


1.407 


1.23 


1.25 


1.27 


2.2 


3.7 


61 


CD2 


1.304 


1.335 


1.388 


2.12 


2.16 


2.22 


3.9 


6.7 


108 


CD 3 


1.313 


1.345 


1.406 


2.79 


2.84 


2.93 


5.1 


8.7 


140 


CD4 


1.555 


1.585 


1.636 


4.54 


4.60 


4.71 


6.5 


11.2 


188 


CD 5 


1.661 


1.682 


1.718 


5.09 


5.14 


5.23 


6.6 


11.6 


197 


CD6 


1.350 


1.383 


1.449 


0.82 


8.36 


0.86 


1.5 


2.5 


40 


CD7 


2.073 


2.078 


2.086 


13.18 


13.21 


13.24 


12.3 


22.0 


389 


CD8 


1.722 


1.780 


1.862 


3.88 


3.99 


4.13 


4.8 


8.3 


141 
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Figure 24. Location of the GW signals from all of our GR models in the h c —fc plane shown relative to the sensitivity curves (giving 
the rms strain noise h IulB ) for the current LIGO detector (left panel) and VIRGO detector (centre panel) and for the advanced LIGO 
detector (right panel). The sources are taken to be at a distance of 10 kpc. Note that some of the models belong to more than one 
sequence. Models for which mode resonance boosts the GW emission are additionally marked with circles. 



(right panel), all for a distance to the source of 10 kpc. As 
a general feature, we note that the upper parts of these 
sensitivity diagrams are occupied by those models whose 
GW signal strength is enhanced by mode resonance. When 
mode resonance is not important, the largest characteristic 
amplitudes occur for those models which are most rapidly 
rotating, unless strong post-bounce damping is in action as 
in model CA1. 

Both in Tableland in Fig. [24] the GW characteristics 
have been evaluated for a total emission time of t{ = 50 ms. 
In Appendix[B]we describe how to obtain the characteristic 
GW strain h c for an arbitrary emission time using simple 
scaling laws. 



6 CONCLUSION 

In this paper, we have described numerical simulations of 
the phase-transition-induced collapse of rotating neutron 
stars to become hybrid quark stars. The simulations were 
made using a code which solves the general relativistic hy- 
drodynamic equations in axisymmetry and within the con- 
formally flat approximation. The initial stellar models were 
taken as being rapidly rotating 7 = 2 polytropes, while dur- 
ing the evolution we used an equation of state composed of 
three parts, depending on density: a normal hadronic-matter 
phase, a mixed phase of normal matter together with de- 
confined quark matter, and a pure quark-matter phase. The 
hadronic component of the equation of state was described 
with an ideal-gas type of equation of state, while for the de- 
confined quark matter phase we used the MIT bag model. 

To validate our code, we first repeated seve ral of the 
Newto nian simulations performed previously by I Lin et all 
|2006t ). We found that the differential rotation which devel- 
ops in these models during the post-bounce phase is almost 
entirely due to strong transient convection which arises be- 
cause the way in which they treated the onset of the phase 



transition leads to an artificial negative specific-entropy gra- 
dient. We argue that their conclusion about there being a 
causal link between the large amplitude post-bounce pulsa- 
tions and the creation of differential rotation was a misin- 
terpretation of the results. We also suggest that a significant 
part of the damping of the pulsations which they observed 
was a consequence of the numerical dissipation present in 
their calculations, rather than being related to differential 
rotation, although a part of this damping was related to the 
other physical processes which we have discussed here and 
which were present in both calculations. 

Having clarified the dynamics of the collapse in the 
Newtonian framework, we then investigated the correspond- 
ing situation within general relativity, using a modified pre- 
scription for triggering the phase transition and initiating 
the collapse. Here we change the equation of state only in 
the regions where the phase transition has taken place and 
leave it unchanged elsewhere. We recognize that our proce- 
dure for describing the phase transition remains extremely 
idealized but we believe that it represents a step forward. 

Despite this difference in the way in which the phase 
transition is triggered, we have not found any major qualita- 
tive differences in the waveforms produced when comparing 
the relativistic simulations with the earlier Newtonian ones. 
Also in the general relativistic simulations, the waveforms 
produced are mainly composed of the fundamental I — 
quasi-radial F-mode and the fundamental / = 2 quadrupo- 
lar 2 /-mode. However, in contrast to the Newtonian models, 
the F-mode is at a lower frequency than the 2 /-mode as 
a consequence of the different density profiles. In addition 
to these modes, a nonlinear self-coupling of the F-mode at 
twice the original frequency, the 2 • F-mode, is strongly ex- 
cited due to the violent nature of the collapse. Although 
qualitatively similar to their Newtonian counterparts, the 
relativistic models exhibit quantitative differences in their 
dynamics. In order to investigate these, we have considered 
a set of 23 different models organized in several sequences. 
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In each of these sequences, only one of the characteristic 
parameters of the models was allowed to vary. 

The main trends observed were as follows. For the se- 
quence with constant initial central rest-mass density, the 
maximum gravitational wave strain |/i| max increases mono- 
tonically with the rotation rate (except for some models 
where the waveform is strongly altered by mode resonances) . 
For the constant rest mass sequence, on the other hand, we 
observe first an increase of |7l| max with the rotation rate and 
then a decrease for very rapid rotation. For the sequence 
with constant rotation period but varying rest mass, we see 
|ft|max increasing monotonically with the rest mass, which is 
a different behaviour from that seen for the Newtonian mod- 
els, where |/i| m ax first increases and then decreases again as 
the rest mass is increased. The reason for this difference may 
simply be that having a decreasing part of the curve would 
require progenitor neutron stars with rest masses beyond the 
upper limit for general-relativistic models. Finally, changing 
the equation of state in the mixed phase has straightforward 
consequences: the central rest-mass density at bounce, the 
amplitude of the post-bounce oscillations, and the maximum 
gravitational wave strain all increase as the overall pressure 
in the mixed phase is reduced. 

Other points to arise include the following. Firstly, 
the influence of rotation on the frequencies of the F 
and 2 f-modes agrees well wit h wha t was found by 
iDimmelmeier. Stergioulas fc Font] (|2006h for pulsations of 
equilibrium neutron stars, suggesting that studies of lin- 
ear pulsations of equilibrium models can (at least in some 
cases) correctly predict the properties of the eigenfrequen- 
cies also when the pulsations are excited in a dynamical 
situation. Secondly, in some models the emission of gravi- 
tational radiation is considerably enhanced by the nonlin- 
ear resonance between the 2 • F-mo&e and the 2 /-mode. 
When the frequencies of these two modes are sufficiently 
close, the weakly radiating quasi-radial 2 • F-mode trans- 
fers some of its kinetic energy to the strongly radiating 
quadrupolar 2 /-mode, leading to a considerable increase 
in |/i| max . Thirdly, we have proposed a simple explana- 
tion for the strong damping of the post-bounce pulsations 
seen for a sub-set of our models. Our analysis reveals that 
that these models are all both rotating close to the Ke- 
pler limit and also undergoing large-amplitude post-bounce 
pulsations, resulting in significant mass shedding from the 
ste llar surface close to the equator . As a lready discussed 
by IStergioulas. Apostolatos fc Fontl (|2004T ). this ejection of 
loosely bound matter is very efficient in the damping quasi- 
radial pulsations. 

In order to assess the prospects for the detection of 
phase-transition-induced collapse events by gravitational 
wave interferometers, special attention has been paid to 
making an accurate calculation of the gravitational wave 
emission resulting from this scenario. We find that the di- 
mensionless gravitational wave strain h from a source at a 
distance of 10 Mpc ranges between 0.1 to 2.4 x 10 -23 for 
all of the models considered and that the total energy emit- 
ted in gravitational waves during first 50 ms of evolution 
is between 10 -6 and 1O _4 M0C 2 , corresponding to 2 x 10 48 
and 2 x 10 50 erg, respectively. These n umbers ar e consi der- 
ably smaller than those presented by I Lin et all (|2006h for 
their Newtonian calculations and so are disappointing for 
the prospects of detecting these sources. The damping times 



for the post-bounce oscillations, as computed from the grav- 
itational radiation waveform, range from 8 to 687 ms for the 
F-mode, and from 18 to 130 ms for the 2 /-mode. For all of 
the models considered, we have also calculated the charac- 
teristic frequency / c , the characteristic strain h c , and the 
signal-to-noise ratio for current and future detectors. For 
current detectors such as LIGO or VIRGO, all of the models 
(except one) have a signal-to-noise ratio above 1 for a source 
at 10 kpc. For the advanced LIGO detector, the signal-to- 
noise ratio rises to well above 10 for a source at 10 kpc for all 
of the models and reaches almost 400 when mode resonance 
is active. However, for detecting sources within the Virgo 
cluster at a distance of 20 Mpc, which is probably neces- 
sary in order to reach a realistic event rate, third generation 
detectors would be needed. 

In conclusion, we note that while our study represents 
an improvement over previous ones, it still lacks a number 
of very important aspects which would be necessary for a 
properly realistic modelling of these objects. Firstly, there 
is the treatment of the phase transformation itself which 
remains extremely crude, containing no detailed picture of 
the transformation of the material, the local heat input or 
the emission of neutrinos or photons. A consistent treatment 
of radiative transfer is likely to be essential for following 
the cooling phase of the newly formed hybrid quark star 
and could highlight that the radiative losses would produce 
differences in the specific-entropy stratification and hence 
trigger real convection. Also, in our discussion, we have been 
considering the phase transition as occurring by means of a 
detonation; the conclusions would be drastically altered if it 
takes place instead via a slow deflagration. A second aspect 
concerns the treatment of the standard neutron-star matter: 
while using a gamma-law equation of state to model this 
can be reasonable for some simplified calculations, it gives 
an extremely poor approximation to the complex physics 
actually occurring in real neutron stars. Thirdly, we have not 
been considering the influence of magnetic fields which are 
not only expected to affect the dynamics, but could also lead 
to a d ramatic modification of the phase transition process 
itself (|Lugones et al.ll2002T ). It is clear that future studies 
will need to take these aspects into account. 
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APPENDIX A: FINITE TIME FOR THE 
INITIAL PHASE-TRANSFORMATION 

In the study by LCCS, it was assumed that the timescale 
of the phase transition is much smaller than the hydrody- 
namic timescale for a NS, which is roughly 0.1 ms. They 
then ignored the finite velocity of the conversion process and 
instead induced the phase transition by instantaneously re- 
placing the initial polytropic EoS by one including the quark 
matter, which gives a lower pressure. 

However, since the timescale for the phase transition 
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even to two-flavour quark matter can be as long as 0.05 ms 
within our picture (which is comparable to the hydrody- 
namic timescale), it is far from clear that treating it as in- 
stantaneous will not significantly affect the subsequent dy- 
namics of the forming HQS. To check on this, we performed 
simulations for a set of representative collapse models in 
which we take the initial phase transformation to occur grad- 
ually over a finite timescale with the values of p nm and p qm 
in Eqs. (|15l I16|l depending on the time t according to 



Phm(t) 



Pqm(t) 



Pc " (Pc - Phm) for t < T CO nv, 

7"conv 

for t ^ Tconv, 



Pi 



Pqm + Pc — Phm 

(Pc — Phm) 

Tconv 

. Pqm 



for t < T CO nv, 



for t > T„ 



(Al) 




Figure Al. Time evolution of the central rest-mass density p c 
for the GR collapse model CA5 with T CO nv = ms (solid line) 
and Tconv = 0.05 ms (dashed line). 



where r conv is the timescale of the initial phase transforma- 
tion and p c is the central rest-mass density. This means that 
Phm(i) starts initially at p c and then evolves linearly with 
time to reach its final value p nm at t = while p qm (£) 

has a similar behaviour shifted by p qm — Phm- For this test 
we set the timescale equal to the approximate upper limit 
for a detonation type phase transition, r coll¥ = 0.05 ms, and 
compared the results obtained with those for the regular 
models (corresponding to T conv — ms). 

We expected that having the initial phase transforma- 
tion occuring gradually in this way would produce a time 
lag of the same order as r conv in the initial contraction and 
give a less violent bounce occuring at a lower density with 
post-bounce pulsations of smaller amplitude than before. All 
of this was indeed the case for the models which we inves- 
tigated, as can be seen in Fig. IA11 where we show the time 
evolution of the central rest-mass density p c for the represen- 
tative model CA5 with T conv = and 0.05 ms. The waveform 
of the emitted GWs remains essentially unaltered except for 
a small reduction in the first large amplitude peaks and the 
expected phase shift (see Fig. IA2[) . As the final HQS is less 
compact in the case of a noninstantaneous phase transition 
due to the less dynamic initial contraction, the frequencies of 
the predominantly excited quasi-normal modes, the F- and 
2 /-mode, are modified only slightly, changing from 1.09 to 
1.12 kHz and from 2.04 to 2.06 kHz, giving relative changes 
of 3% and 1%, respectively. 

The result that the differences seen when the finite 
Tconv is introduced are so small is a direct consequence of 
Tconv = 0.05 ms still being about an order of magnitude 
smaller than the collapse timescale in the case of an in- 
stantaneous initial phase transformation (although it is of 
the same order as the dynamical timescale). The collapse 
timescale can be approximated by the time of bounce, which 
is 0.40 ms for model CA5 and takes similar values for all of 
the other models. Since the choice T conv — 0.05 ms is an 
upper limit within our picture, we conclude that the taking 
Tconv = ms for our regular collapse models was reasonable, 
giving values for quantities such as p c ,b, |ft|max, /f and f2f 
that are upper, but close, limits for a detonation type phase- 
transition- induced collapse of a NS to a HQS. 
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Figure A2. Time evolution of the GW strain h at a distance of 
10 Mpc for the GR collapse model CA5 with r conv = ms (solid 
line) and T con v = 0.05 ms (dashed line). 
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Figure A3. Power spectrum h (in arbitrary units) of the GW 
strain h for the GR collapse model CA5 with T con v = ms (solid 
line) and T CO nv = 0.05 ms (dashed line). 
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APPENDIX B: CALCULATION OF THE 
DAMPING TIMES OF THE WAVEFORMS 

In order to determine the damping times for the GW emis- 
sion we apply a curve fitting method to the numerically- 
obtained waveform: the waveform is fitted by a series of 
damped sinusoids, 

n 

Am =J2 Ai fr ' /T * cos ( 27r /*i + 00- (Bl) 

The parameters of this series (the damping times n, the 
amplitudes Ai and the phases (pi) are fixed by finding the 
best fitting curve. Depending on the model, we use n = 4-5 
terms for the fitting procedure. 

Note that this method of determining the damping 
times is more general than that presented in LCCS, where 
the series consisted of only two terms, so that only the fun- 
damental I — and I — 2 modes are taken into account, 
which leads to some inaccuracies at early post-bounce times 
when higher frequency modes in the GW signal are present. 
Also, they take Ai = A2 and (pi = 0, which may not be 
the case in general. With our approach we are able fit the 
original waveform with the correlation coefficient between 
the numerical data and Eq. (|B1|) exceeding 0.99 for all mod- 
els, whereas the method used by LCCS gives a correlation 
coefficient of less than 0.9 for some models. 

Assuming that the GW signal is essentially a linear com- 
bination of the F-mode and the 2 /-mode, approximating 
each of them as a damped sinusoid and using knowledge of 
the mode frequencies fp and faf, the amplitudes A\ and A2, 
and the phases cj>\ and 4>2, one can obtain the value of the 
characteristic GW strain h c for an arbitrary emission time 
using simple scaling laws. For a single damped sinusoid, 

h = ho e' t/T sin (2irft - <p) (B2) 

where ho is the amplitude, r is the damping timescale, / is 
the frequency and <f> is the phase. If h c is known for some 
emission time t{ , then its value for a multiple n of the original 
emission time can be calculated as 

fec(ntf) = fe c (tf) J 1 _ e _ 2t{/T , (B3) 

provided that / _1 <C r (which is fulfilled for most of our 
models) and the power spectral density Sh of the detector is 
reasonably constant in the vicinity of /. For an undamped 
sinusoid with r = 00, Eq. ()B3|) gives h c (nt{) — y/nh c as 
expected; in other words, h c scales like the square root of 
the number of cycles in the GW signal. In the limit of infinite 
emission time, but with finite r, the exponential damping of 
the signal results in a finite value for the total characteristic 
GW strain, 

h c (t=oo) = h c (t c ) 1 (B4) 

VI — e~ 2t f/ T 

For many of our models, we have f-mode damping 
times tf which are much longer than t{ . On these timescales, 
other damping mechanisms such as physical viscosity or 
GW back-reaction (which are not taken into account in our 
study) could become important. We therefore do not give 
values of the total characteristic GW strain h c (t = 00) for 
our models. 



